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Summary 


In August 2014, the Rosetta spacecraft completed its ten year journey when it arrived 
at its target destination, the comet 67P/Churyumov-Gerasimenko. The Rosetta mission 
was a flagship endeavour for the European Space Agency as it was the first time that any 
spacecraft had rendezvoused with a small solar system body for long period of time (two 
years) and also the first time that a lander had been placed onto the surface of a comet. 
In September 2016, the mission came to an end when Rosetta descended onto the surface 
for one final close up look at the surface. 

This thesis uses data from one of the instruments on this ground-breaking mission: 
the Microwave Instrument for the Rosetta Orbiter (MIRO). MIRO was a spectrometer 
operating at millimetre and submillimetre frequencies and enabled the detection of several 
volatile species including water. Using the spectroscopic observations of the water lines, 
I investigated the behaviour of comet 67P/Churyumov-Gerasimenko relating to the gas 
coma, mass loss, spatial outgassing and the nucleus surface. Since comets are thought 
to be pristine building blocks left over from the formation of the solar system, we hope 
that by studying them, we can learn about the conditions from which other solar system 
bodies originate. 

Firstly, I used the line areas of the H}°O and H3°O spectral lines to measure the change 
in the local water production rate from August 2014 to April 2016. Lookup tables made 
from a Haser model show how the measured Doppler shift velocity, the continuum tem- 
perature and the line area ratio can give the column density for each observation and 
thus the water production rate. The maximum production rate calculated from the MIRO 
observations was (1.42+0.51)x108 molec/s, found on August 29, 2015, and integrating 
under all the data points gave a total water mass loss of (2.4 + 1.1)x10° kg for this ap- 
parition. By making assumptions about the dust-to-gas ratio and the comet mass, the 
total mass loss was estimated to be (1.2 + 0.6)x10!° kg, or 0.12 + 0.06 %. The spatial 
resolution of MIRO allowed for each measurement to be assigned to a region on the nu- 
cleus. The regions on the southern hemisphere appeared to be the origins of the highest 
production rates, in particular the regions Neith, Wosret and Bes. Finally, the data show 
that the production rate peak is offset by 22-46 days after perihelion and that the pre- and 
post-perihelion slopes followed power laws of Q c r,?5*9? and Q « r,*2*9?. respectively. 

Following this, I performed numerical modelling to investigate how nucleus shape, 
spin axis orientation and activity distribution affect the water production rate curves. I 
found that it is impossible to disentangle these effects from each other by only looking 
at the change in the production rate and that the pre- and post-perihelion slopes are also 
functions of heliocentric distance. It is therefore difficult to derive quantitative constraints 
on surface area ice fraction and active regions from the water production rate curve unless 
the shape, orientation and activity of the nucleus are well established. In addition, it may 
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not be meaningful to compare the water production rate curves of different comets at 
different heliocentric distances. 

I used the measured continuum temperatures from MIRO to derive properties of the 
nucleus in the third part of this work. I utilised an insolation driven thermal model to 
derive the temperature gradient in the upper layers of the comet surface and a radiative 
transfer model to reproduce the MIRO continuum measurements. In conclusion, a low 
value was derived for the thermal inertia in the surface layers of 67P with an upper bound 
estimated to be 80*7 JK^! m ?s 9? for most of the MIRO measurements. A low value for 
the average thermal inertia over the entire surface would be consistent with the majority 
of reported calculated values for 67P. 

In the future, the retrieval of coma parameters from the MIRO spectra will become 
an important avenue of investigation. Using inverse methods, the behaviour of the gas 
temperature, expansion velocity and molecular number density profiles can be extracted 
from the spectral lines. This will be important for assessing and characterising the ac- 
tivity around the nucleus which is observed but not so well understood. In addition, we 
may learn more about the physics of the coma from the inversion solutions, such as the 
presence of a Knudsen layer, photolytic heating or extended sources. At the end of my 
thesis, the application of inverse methods to the MIRO data is described. 

It is still uncertain to what extent 67P is pristine and primordial but the results of my 
thesis imply that the surface layer over a couple of metres must be processed by the sun 
during perihelion passage. In this region, icy volatiles sublimate as the thermal wave 
propagates into the surface, expelling dust and refractory material away from the nucleus 
along with the gas. In the most productive regions on the south, the erosion is expected 
to reach down to a few metres in depth. The proposed CEASAR sample return mission 
would revisit 67P and the work for my thesis suggests that the drill must go beyond the 
surface layer, perhaps as much as 10 times the diurnal skin depth, in order to reach any 
potentially pristine material. 
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Zusammenfassung 


Im August 2014 hat die Raumsonde Rosetta ihre zehnjährige Reise beendet, als sie am 
Ziel, dem Kometen 67P/Churyumov-Gerasimenko angekommen ist. Die Rosetta-Mission 
war eine Flaggschiff-Mission der europäischen Weltraumagentur, da die Raumsonde zum 
ersten Mal mit einem kleinem Körper im Sonnensystem für einen längeren Zeitraum 
(zwei Jahre) in Kontakt war und eine Landeeinheit zum ersten Mal auf einer Kome- 
tenoberfläche platziert wurde. Als Rosetta im September 2016 für einen letzten detail- 
lierten Blick zur Kometenoberfläche herunter gefahren ist, fand die Mission ihr Ende. 

Diese Doktorarbeit verwendet Daten eines der Instrumente dieser bahnbrechenden 
Mission: das Mikrowellen-Instrument für den Rosetta Orbiter (MIRO). MIRO war ein 
Spektrometer, welches im Millimeter- und Submillimeter-Frequenzbereich arbeitete und 
die Entdeckung mehreren flüchtigen Spezies einschließlich Wasser ermöglicht hat. Bei 
der Anwendung der spektroskopischen Beobachtungen der Wasserlinien habe ich das Ver- 
halten des Kometens 67P/Churyumov-Gerasimenko in Bezug auf die Koma, den Massen- 
verlust, räumliche Ausgasungen und die Kernoberfläche untersucht. Da Kometen als vor- 
malige Bausteine aus der Zeit der Entstehung des Sonnensystems angesehen werden, er- 
hoffen wir, dass wir bei deren Untersuchung über die Entstehungsbedingungen anderer 
Körper im Sonnensystem lernen können. 

Zuerst habe ich die Linienintegrale der Spektralinien H160 zur Messung der Verän- 
derung der lokalen Wasserproduktionsrate zwischen August 2014 und August 2016 unter- 
sucht. Mit Hilfe von Look-Up Tabellen aus dem Haser Model kann von den gemessenen 
Doppler-Geschwindigkeiten, der Kontinuums-Temperatur und dem Verhältnis der Lin- 
ienintegrale auf die Säulendichte und somit die Wasserproduktionsrate geschlossen wer- 
den. Die von den MIRO-Beobachtungen gerechnete höchste Wasserproduktionsrate war 
(1.42+0.51)x10°® molec/s, die am 29. August 2015 gefunden wurde. Beim Intergri- 
eren aller Datenpunkte hat sich für dieses Ereignis ein Gesamtwasserverlust von (2.4 + 
1.1)x10? kg ergeben. Mit Hilfe von Annahmen über das Staub-Gas-Verhältnis und der 
Kometenmasse, wird der Gesamtverlsut auf (1.2 + 0.6)x10!° kg oder 0.12 + 0.06 % der 
Kometenmasse geschätzt. Dabei hat es die räumliche Auflösung von MIRO ermöglicht, 
dass jeder Messung eine Region am Kometenkern zugewiesen werden konnte. Die Regio- 
nen an der Südhalbkügel, insbesondere die Regionen Neith, Wosret und Bes erschienen 
als die Ursprünge der höchsten Produktionsraten. Letztendlich zeigen die Daten, dass 
die maximale Produktionsrate um 22 bis 46 Tage zum Perihelion versetzt ist und die 
Produktionsraten vor und nach dem Perihelion den Potenzgesetzen Q « r;?9*9? bzw. 
Q « r,^2*9? folgen. 

Im Anschluss daran habe ich numerische Modellierungen durchgeführt, um zu un- 
tersuchen, wie die Kometenform, die Ausrichtung der Rotationsachse und die Aktiv- 
itätsverteilung die Wasserproduktionsrate prägen. Ich habe festgestellt, dass nicht möglich 
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ist diese Effekte getrennt voneinander zu betrachten, wenn lediglich die Veränderung der 
Produktionsrate betrachtet wird. Zudem sind die Produktionsraten vor und nach Perihe- 
lion abhängig von der heliozentrischen Distanz. Deswegen ist es schwierig, die quan- 
titativen Beschränkungen auf der Oberfläche der eisbedeckten Gebiete und der aktiven 
Regionen von der Wasserproduktionsratekurve herzuleiten, außer wenn die Kernform, 
die Kernausrichtung und die Kernaktivität gut bestimmt sind. Zudem könnte es nicht sin- 
nvoll sein, die Wasserproduktionsratenkurven verschiedender Kometen bei verschiedenen 
heliozentrischen Distanzen zu vergleichen. 

Im dritten Teil dieser Arbeit habe ich die gemessene Kontinuumstemperatur von MIRO 
zur Herleitung der Kerneigenschaften verwendet. Ich habe ein von Einstrahlung reg- 
uliertes thermisches Modell zur Derivation das Temperatursgradienten in den oberen 
Schichten der Kometenoberfläche verwendet, sowie ein Strahlungstransportmodell zur 
Wiederherstellung der MIRO-Kontinuumsmessungen. Damit wurde ein niedriger Wert 
für die thermische Trägheit in den Oberflächenschichten von 67P gefunden. Der Wert 
hatte eine geschätzte obere Grenze von 80*%) JK !m ?s 9? für die Mehrheit der MIRO- 
Messungen. Ein niedriger Wert für die durchschnittliche thermische Trägheit über die 
Gesamte Oberfläche, würde konsistent mit dem Großteil der für 67P errechneten Werte 
sein. 

In der Zukunft wird die Gewinnung von Kometenparametern von den MIRO-Spektren 
eine wichtige Untersuchungsrichtung. Bei der Verwendung von inversen Methoden kön- 
nen das Verhalten der Gastemperatur, der Ausbreitungsgeschwindigkeit und der 
Molekülanzahldichteprofile von den Spektrallinien extrahiert werden. Dieses Ergebnis 
wird wichtig für die Beurteilung und Charakterisierung der Aktivität um den Kern sein, 
die beobachtet wird , aber noch nicht gut verstanden ist. Zusätzlich können wir durch In- 
versionslósungen unser Wissen der Physik der Koma erweitern, wie die Gegenwart einer 
Knudsen-Schicht, photolytische Erwärmung oder erweiterte Quellen. Am Ende meiner 
Doktorarbeit wird die Anwendung inverser Methoden auf die MIRO-Daten beschrieben. 

Es ist immer noch unklar, inwieweit 67P unberührt und unberührt ist, aber die Ergeb- 
nisse meiner Arbeit deuten an, dass die Oberflächenschicht über eine Dicke von weni- 
gen Metern während des Periheldurchgangs durch die Sonne bearbeitet wurde. In dieser 
Region sublimieren Eisvolatile während die thermische Welle sich in die Oberflächen- 
schichten ausbreitet, wobei sie Staub und hitzebständiges Material, sowie das Gas vom 
Kern ausgestoßen werden. In den hochaktiven Region im Süden ist die Erosion vo- 
raussichtlich ein paar Meter tief. Die vorgeschlagene CEASAR-Mission, die eine Probe 
liefern soll, würde 67P wieder besuchen. Meine Arbeit weist darauf hin, dass das Bohrg- 
erät zur Erreichung eventuell unberührtes Materials durch die Oberflächenschicht weiter 
hinaus, vielleicht so tief wie zehnmal der tageszyklischen Eindringtiefe, bohren muss. 
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1.1 Comets 


Comets are small solar system bodies made from icy volatile and rocky refractory mate- 
rial. They are typically a few kilometres in size (1-10 km), with the Great Comet of 1577 
(C/1577 V1) and Comet Sarabat of 1729 (C/1729 P1) thought to be the largest with diam- 
eters from 100-300 km (Clube et al. 1996), whereas the bi-lobed comet 103P/Hartley has 
a maximum length measuring only 2.33 km (A’ Hearn et al. 2011). Due to their small size, 
comets are usually invisible when they are far away from the Sun and appear as only point 
sources, indistinguishable from the background stars. However, when they come close to 
the sun, comets become much brighter and develop observable comae and tails, as shown 
in Figure 1.1. The volatile content within the nucleus sublimates as a result of the increase 
in received solar flux as the heliocentric distance decreases and the sublimated gas forms a 
coma, a tenuous atmosphere around the comet. The coma also contains dust grains which 
are lifted from the nucleus by the outflowing gas, although this historical idea seems to 
be somewhat more complicated as the exact mechanism of lifting dust is unknown - gas 
drag appears to be insufficient to break the cohesion force between grains on the surface 
(Skorov et al. 2017). The size of the coma depends upon how it is defined (Rodgers et al. 
2004): the collisional region may be 10-100 km in size and the H2O coma can be 10? km 
in radius, whilst the atomic hydrogen coma can measure as much as 10’ km. However, 
the material in the coma is not gravitationally bound due to the small size of comets. An 
ion tail grows as ionised gas and fine grained dust is pushed away from the nucleus, and 
a dust tail forms from solar radiation pressure. Both tails always point away from the 
sun, but the ion tail is always in a radial direction from the sun due to interactions with 
the solar wind, and the dust tail is curved, as the small dust grains start to move on their 
own Keplerian orbits (Karttunen et al. 2016). The tails and the coma become optically 
observable due to the dust material within them reflecting light. 

Comets are almost ubiquitous throughout the solar system, with some comets having 
orbits which keep them close to the terrestrial and giant planets, others on long period 
orbits beyond Neptune (>30 AU) which take them through the scattered disk and the 
Kuiper Belt, and some located in the furthest reaches of the Oort Cloud (>1000 AU), the 
spherical shell that encompasses the solar system. However, the most distant observation 
of a comet was of Hale-Bopp at 32 AU (Szabo et al. 2012) and due to their dimness, 
it is difficult to observe any comets beyond that in the scattered disk or the Oort Cloud. 
The evidence for comets in these regions comes from the reconstructed orbits of comets 
which enter the inner solar system. The most distant comets are the nearly isotropic 
comets (NICs) which have semi-major axes from ~40 AU to over 10,000 AU and they 
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Figure 1.1: C/1911 Ol (Brooks), discovered by William Robert Brooks. (Credit: Paul 
Anderson, 1911) The image shows the large comet tail forming behind the nucleus. 


reside beyond Neptune, in the scattered disk, the Kuiper Belt and the Oort Cloud (Levison 
1996, Levison and Dones 2014). There are three types of NICs, new comets, which 
are the furthest away, returning comets, and Halley-type comets. Ecliptic comets are 
found in the orbital plane and have semi-major axes less than 40 AU, bringing them 
under the gravitational influence of the planets. One particular class of ecliptic comets are 
called Jupiter-Family Comets. These comets cross the orbit of Jupiter and feel its strong 
gravitational attraction (Duncan et al. 2004, Dones et al. 2004). 

The dynamical lifetimes of comets are short compared to the age of the solar system 
and so they must have been put on their current orbits relatively recently and come from 
from some other reservoir (Duncan et al. 2004). Two pieces of observational evidence 
suggest that NICs and the ecliptic comets originate from different reservoirs (Levison and 
Dones 2014). The first is that while some comets change minor class (i.e. Jupiter-Family 
comet to Chiron-type, or new to returning comet), comets do not generally change major 
class (NIC to ecliptic) (Levison 1996). In addition, the inclination distribution for NICs 
is very isotropic, suggesting they come from all directions (hence the name for this class) 
(Dones et al. 2004), but the inclination distribution for ecliptic comets is much narrower 
suggesting they come from a different region (Levison and Dones 2014). Even Halley- 
type comets have a more confined range of inclinations than new and returning comets. 
It is generally thought then that NICs come from from the Oort Cloud at the edge of the 
solar system, coming into the inner solar system after gravitational perturbations from 
outside influences; and that ecliptic comets have origins in the scattered disk, a region 


14 


1.1 Comets 





that extends from Neptune to the inner edge of the Oort Cloud, where numerous trans- 
Neptunian objects have been observed. Halley-type comets appear to originate from an 
intermediate region where the spherical Oort cloud is becoming more disk-like. 

However, it is also thought that comets did not initially form in these reservoirs, and 
that they have a complex dynamical history due to the supposed migration of the planets 
during the early formation of the solar system. It is thought that most comets should have 
formed around and beyond the snow line where the giant planets are currently located 
(Dones et al. 2004, 2015). In the early evolution of the solar system, as the swirling disk 
of gas rotated around the proto-Sun, small planetesimals started to form and grow in size. 
Some of these objects coalesced to form the planets whereas others where slowly kicked 
out of the inner solar system by planetary interactions, increasing their semi major axis 
and perihelion distance. These planetesimals which were kicked out into the scattered 
disk and the Oort Cloud are the comets which we can now see today. 

Through spectroscopic observations and satellite flybys, the composition of comet 
nuclei is slowly being revealed. In 1950, a revolutionary time for cometary science, it 
was first proposed that comets are dirty snowballs (Whipple 1950) made from volatile 
ices and dusty refractory material. Our knowledge of the volatile component arises mainly 
from spectroscopy, which has allowed dozens of molecules to be identified in the coma 
of comets (Bockelée-Morvan 2011). The main component is water, accounting for up to 
80% of the volatile content. The abundance of other species are usually given as relative 
to the water content. After water, CO and CO, are two of the more abundant molecules 
observed in cometary comae but the abundances in the coma vary from comet to comet. 
For example, in most comets so far observed, CO»; dominates CO, however for comet 
C/2006 W3 (Christensen), the situation appears to be reversed and CO dominates CO, 
(Ootsubo et al. 2012). Deriving classification based on the ratios of volatile abundances 
is very difficult though due to the high variability in the abundances, for example, the 
CO;/H50 ratio can vary with heliocentric distance and with latitude on the surface (Mall 
et al. 2016). There is also further evidence that the molecular abundances in the coma 
are not constant as seen with comet 67P/Churyumov-Gerasimenko, where the number 
density ratio of CO;/H5O has been observed to vary between 2.5926 and 80% (Fougere 
et al. 2016) as the heliocentric distance changes and as observations are made on different 
parts of the nucleus. Abundances can be measured for a single observation but this is only 
a single snapshot of the overall complicated outgassing behaviour. Volatile molecules 
have different sublimation temperatures and so the outgassing of certain ices switches 
on at different heliocentric distances. All of this must be considered when interpreting 
volatile abundance ratios. 

Other molecules observed in relatively high abundances include CH3OH, CH,4, H2S 
and NH3. These species, in addition to CO and COs, are usually observed to have abun- 
dances of 0.01% - 20% relative to water. Figure 1.2 lists a number of observed volatile 
species and their relative abundances to water. 

There are ongoing questions as to whether the spectroscopic observations made of 
the coma are representative of the nucleus, and even how much the nucleus may have 
evolved since its formation. Comets may experience some evolution far away from the 
sun in the Oort Cloud as well as during close perihelion passages when the volatile mate- 
rial sublimates (Bockelée-Morvan 2011). Evidence from 73P/Schwassmann-Wachmann 
3, which started to disintegrate as it entered the inner solar system, shows that the dif- 
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Figure 1.2: Abundance ratios relative to water for a variety of molecules observed in 
cometary comae (based on Fig. 1 of (Bockelée-Morvan 2011)). 


ferent cometary fragments have similar chemical abundances, suggesting that chemical 
abundances could be primordial signatures of formation rather than from evolutionary 
processes over the comets lifetime (Cochran et al. 2015). 

Information on the refractory part of comets mainly arises from cometary space mis- 
sions (Bockelée-Morvan 2011). The dust collected from around comet 1P/Halley was 
analysed using mass spectrometry to determine the composition. Elements such as Mg, 
Si, Fe and Ca were found, as well as organic material including CHON grains. 

The study of comets is thought to be very important as they have profound impli- 
cations for our solar system. Far from the sun where the temperatures are low, comets 
should remain relatively unprocessed over their lifetimes and the interior material must 
be relatively pristine and primordial. If comets do originate from the very early formation 
of the solar system and were ejected to the furthest reaches of the sun, then they may hold 
clues to the conditions and evolution history of our solar system and could shed light on 
what was happening approximately 4.5 billions of years ago. 


1.2 Observational history: from the stone age to the space 
age 


Remarkably, it is argued that the earliest recorded sightings of comets come from the 
latter half of the Shang (or Yin) Dynasty which ruled in Ancient China between 1500 - 
1050 BC. The observations were recorded on animal bones and turtle shells and they are 
known as the Oracle Bones, containing records of solar eclipses, lunar eclipses, planets 
and comets (Zhen-Tao et al. 1995). A few hundred years later, the Babylonians also 
recorded observations of comets in their Astronomical Diaries, which covered the period 
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from 750 BC - 75 AD, and wondered about their origins, speculating that these remarkable 
objects in the sky were obscure planets (Stephenson et al. 1985). The Babylonian records 
also include mentions of probable observations of comet 1P/Halley in 164 BC and 87 BC. 
Ephorus, the Greek historian, believed that comets were the result of collisions between 
stars after observing the Great Comet of 371 BC, which is now suspected to have been a 
Sun-grazing comet that broke apart during this apparition (Seargent 1982). 

However, without a proper explanation, other people believed that, rather than be- 
ing astronomical bodies, comets were actually meteorological phenomenon (Seargent 
1982). Numerous famous figures throughout history, including Pliny the Elder, Ptolemy 
and Galileo, thought that comets were atmospheric manifestations like reflections of the 
Sun or features from volcanic eruptions. 

In addition, comets were also mistakenly regarded as portents of doom by some peo- 
ple, heralding a natural disaster, the death of a monarch or attacks by heavenly forces. 
Across Europe, the Middle East and Asia, the appearance of a comet was treated with 
suspicion and superstition. The 164 BC apparition of 1P/Halley coincided with the death 
of Antiochus IV, the ruler of the Selucid Empire, as well as the expulsion of the Ptole- 
maic ruler Philometor from his throne (incidentally, both of their successors were replaced 
within two years) (Wolters 1993). In 44 BC, the apparition of a comet coincided with the 
death of Julius Caesar, as written about by William Shakespeare (Seargent 1982). The 
Bayueux Tapestry, a medieval piece of artwork depicting the successful conquer of Eng- 
land by William the King of Normandy, also envisions a comet in the background, an 
omen of the end of the reign of the English King Harald. Giotto di Bondone, the 14" 
century painter, was so taken by the mystical elements of comets that he chose to include 
a comet in his painting, The Adoration of the Magi, showing the birth of Jesus Christ in a 
nativity scene (Crovisier et al. 2000). 

It was not until Tycho Brahe made astronomical observations of the solar system dur- 
ing the 16* Century that the idea of comets as astronomical bodies returned to prominence 
(Seargent 1982). Brahe made meticulous measurements of the Great Comet of 1577 from 
his observatory in Hven, Sweden, in conjunction with Thaddeus Hagecius in Prague, and 
determined the parallax for the comet (Christianson and Brahe 1979). In their observa- 
tions, the comet was approximately in the same position on the sky for both of them but 
the Moon was not, leading Brahe to conclude that the comet must be further away from 
Earth than the Moon and therefore not an atmospheric phenomenon. A student of Brahe, 
Johannes Kepler, confirmed this finding by measuring the parallax of the comets of 1607 
and 1618 (Seargent 1982). 

The exact nature of comets was still unknown in this time but Edmond Halley shed 
some more light on the nature of cometary objects when he observed comets in 1680 
and 1682 (Seargent 1982). After calculating the orbit in collaboration with Isaac Newton 
(Gribbin 2004), he speculated that the second one may have been the same comet which 
had appeared in 1607 and 1531, and that comets had a periodic nature. He predicted that 
the comet would return in 1758. When the comet reappeared, Halley's prediction proved 
to be correct, although he did not live to see this as he died in 1742, sixteen years before 
the apparition. In addition, Halley thought, much like Kepler before him, that the comets 
which they had seen were only one of many in existence. 

It was not until the year 1950 that cometary science took a huge leap forward with the 
publication of three seminal papers which underpin our basic understanding of comets 
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(Festou et al. 2004). The first was the "dirty snowball" view of comets proposed by Fred 
Whipple (Whipple 1950), which postulated that comets were icy conglomerates contain- 
ing a mixture of volatile material that sublimated in increasing quantities as the comet 
approached the sun and releasing dust grains from the surface. It was an important idea 
which encompassed many of the observed aspects of comets which had been documented 
in earlier work (Festou et al. 2004), such as: the presence of molecular daughter species 
(CO*, CN) from more stable ices (CO, CO;); the large production rates of comets; jet-like 
structures; the influence of nongravitational forces on their motion; and meteor streams. 
The work of Whipple (1950) was widely accepted at the time and it is still of major 
importance today. 

The second novel idea was Jan Oort's idea (Oort et al. 1950) of a halo surrounding 
the solar system containing a huge number of comets: the Oort Cloud. This idea evolved 
from the distribution of the semi major axis of comets. Over time our understanding of 
the Oort Cloud has improved but the work of Oort has informed our current understanding 
and helped define the taxonomic classes described in Section 1.1. 

The connection between the comet tail and the solar wind was the third finding in 
this time frame from Biermann (1951). It was already suspected that magnetic activity 
on Earth was influenced by solar flares (Festou et al. 2004), but it was Biermann's work 
which hypothesised that the ions in a cometary tails were interacting with electrically 
charged particles of the solar wind. 

Clearly, our knowledge of comets has improved over the centuries, and we now know 
that comets are not meteorological or prescient objects at all, but bodies in our solar 
system. With the advent of space exploration in the latter half of the 20" century, it 
became possible to fly instruments to distant comets and study them close up for the very 
first time. 

Studying objects deep in the Oort Cloud, Kuiper Belt or the scattered disk is not easy 
though, and going to them for close up exploration is almost impossible. If we want to 
study comets, we must then look a little closer to home. The Jupiter Family Comets, 
although processed by their close approaches to the Sun, may still contain some of their 
pristine nature and by studying them in great detail, we can try to answer one of the 
greatest mysteries in space science: how the solar system came to be. 

The history of cometary spacecraft exploration is only 33 years old and in total, only 
eight different comets have been visited by ten different spacecraft (Schulz 2009). The 
first was the NASA spacecraft, the International Cometary Explorer (ICE), which visited 
comet 21P/Giacobini-Zinner in 1985 and made measurements of particles, waves and 
fields. ICE probed the interaction between the comet and the solar wind, confirming the 
Alfven model structure of a comet tail (two lobes of opposite polarity magnetic field lines 
separated by a current sheet), detecting water and carbon monoxide ions with the water 
group ions dominating, and measuring collisions of dust particles with the spacecraft 
(Von Rosenvinge et al. 1986). 

A year later, a collection of five spacecraft, unofficially known as the Halley Armada, 
visited comet 1P/Halley during its apparition in 1986. These were: the Russian probes, 
VeGa 1 and VeGa 2; Sakigake and Suisei, built by the Institute of Space and Astronautical 
Science (ISAS, now known as the Japanese Aerospance Exploration Agency, JAXA); and 
Giotto from the European Space Agency (ESA) (Festou et al. 2004). VeGa 1 sent back 
the first images of a comet nucleus (Sagdeev et al. 1986) and measured the temperature of 
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Figure 1.3: Spacecraft images of comets. (1.3a) the nucleus of 1P/Halley as seen by 
the ESA Giotto probe in 1986 (Credit: ESA, 2000-2008 © European Space Agency), 
(1.3b) the nucleus of 81P/Wild, taken by the Stardust mission in 2004 (Credit: Courtesy 
NASA/JPL-Caltech.), (1.3c) the nucleus of 9P/Tempel, seen from the Deep Impact mis- 
sion in 2005 (Credit: Courtesy NASA/JPL-Caltech.), (1.3d) the nucleus of 103P/Hartley 
as seen by EPOXI in 2011 (Credit: Courtesy NASA/JPL-Caltech.) 
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the surface to be 300-400 K, much warmer than expected for an icy body (Combes et al. 
1986, Emerich et al. 1988). This implied there must be a layer of dust material at the sur- 
face. VeGa 2 made a closer approach to the surface and enabled the determination of the 
physical size of the nucleus (Sagdeev et al. 1986). Suisei took measurements of 1P/Halley 
with an ultraviolet imager and confirmed the hard surface of the nucleus, as well as de- 
termining some solar wind parameters (Kaneda et al. 1986). Observations from Sakigake 
lead to the proposal of an Interplanetary Magnetic Field-comet interaction model (Saito 
et al. 1986). Finally, Giotto encountered Halley, getting even closer to the nucleus, and 
became the first instrument to make mutlicolour images of a comet (Thomas and Keller 
1988), see Figure 1.3a. This revealed the jet-like structures made from fine dust as well 
as morphological features (Schwarz et al. 1988). Giotto also helped to constrain predom- 
inance of water in the coma (Krankowsky et al. 1986), and the size and composition of 
the ejected dusty material, including the presence of CHON elements (Clark et al. 1988). 

In 1992, the Giotto Extended Mission visited the comet 26P/Grigg-Skjellerup, making 
the closest approach of any comet until the 2014 Rosetta mission (Schulz 2009). Despite 
the failure of the camera, which was damaged and inoperable after its encounter with 
Halley, it was found that comet 26P/Grigg-Skjellerup is much smaller than 1P/Halley and 
relatively inactive in comparison. 

Fifteen years passed before the next cometary mission when Deep Space 1 travelled to 
comet 19P/Borrelly. Although the main goal of the mission was a flyby of asteroid Baille, 
the Deep Space encounter with Borrelly was incredibly successful, making images of a 
Jupiter Family Comet for the first time, and at a higher resolution than ever before (Boice 
et al. 2000). The images revealed an exceptionally dark surface, with albedo variations 
from 0.007 - 0.035, and a hot, dry surface with a sublimating area less than 10% of 
the whole nucleus. The PEPE plasma instrument also observed inhomogeneities in the 
plasma data around the comet (Richter et al. 2011). 

Comet 81P/Wild was the subject of the first ever comet sample return mission, with 
the Stardust spacecraft making a rendezvous with the comet in 2004, before completing 
its mission in January 2006 when the samples reached Earth (Brownlee 2014). Stardust 
contained several deposits of low density silica aerogel to collect dust particles ejected 
from the surface of the comet and stored them for return and analysis in a laboratory. 
The returned particles shed a lot of light on the early history of the Solar System. The 
dust component was similar in material to meteorites and contained a low abundance of 
material from the presolar nebula. This implies that 81P/Wild formed during the very 
early stage of the solar system and from material which condensed at different heliocen- 
tric distances, from the hot interior where the rocky components condensed to the cold 
outer regions where the supervolatile material formed. The images from Stardust showed 
81P/Wild to be a 5 km, round body, with many ancient surface features despite its active 
nature (Brownlee et al. 2004), as shown in Figure 1.3b. 

The next cometary mission was the visit of Deep Impact to 9P/Tempel in 2005. The 
unique aim of this mission was to study the interior composition of a comet by firing 
an impactor, a 360 kg projectile, at the nucleus and thus analyse the resulting exposed 
debris for pristine material (Hampton et al. 2005). The observations from Deep Impact 
revealed several things about the comet (A’ Hearn et al. 2005): the nucleus has a density of 
approximately 600 kg m^? and a diameter of about 3 km; there was an observed increase 
in organic material relative to water during and after the impact; the upper layer of nucleus 
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is made of fine particles (1-100 um) with negligible strength and these particles may be 
present to decimetre scale depths; the surface is in equilibrium with sunlight; and ice 
does not appear to be on the surface so the observed outgassing must come from volatile 
material near the surface. A Stardust image of 9P/Tempel is shown in Figure 1.3c 

Renamed EPOXI, the Deep Impact spacecraft was requisitioned for two further mis- 
sions. The first was to perform a fly-by of another comet, and the second aim was to use 
its imaging cameras as a space observatory to observe exoplanets and comets. In Novem- 
ber 2010, EPOXI made a close fly-by of comet 103P/Hartley, a very active comet despite 
being quite small, depicted in Figure 1.3d. In fact, the comet was producing more water 
than would be possible from an object of its size (A’Hearn et al. 2011). CO, must there- 
fore have been an important driver of activity, releasing icy dust grains from the surface 
to enhance the water production. From the images, it appeared as if most of the water 
outgassing was coming from the waist of the elongated comet. During the observatory 
phase, EPOXI made observations of Garradd (C/2009 P1) (Feaga et al. 2013) and C/2012 
S1 (ISON) (Farnham et al. 2017), measuring the volatile production rate of these comets. 
The mission terminated in 2013 after contact with the spacecraft was lost. 

The Stardust mission was extended (called Stardust-NExT) and made a flyby of 
9P/Tempel in 2011, enabling direct comparisons of observations from Deep Impact (Vev- 
erka et al. 2013). The additional observations from Stardust-NExT enabled the coverage 
of the nucleus to be extended to two-thirds of the surface in combination with the Deep 
Impact images. The surface showed evidence of layering, deep pits in the nucleus, and 
some changes in small regions since the 2005 fly-by, although the majority of the surface 
was relatively unchanged. The jet activity also appeared lower than before. 

The most recent exploration of a comet occurred only three years later with the arrival 
of the Rosetta spacecraft at the comet 67P/Churyumov-Gerasimenko. 


1.3 The Rosetta mission 


The Rosetta mission was first envisaged in 1991 as a sample return mission as part of the 
European Space Agency Horizon 2000 programme (Huber and Schwehm 1991). ESA 
approved the Rosetta mission in 1993. The sample return aspect of the mission was 
later dropped but many of the key aims remained unchanged: to probe the origins of 
cometary formation and evolution, as well as to learn about the origins of the solar system 
(Glassmeier et al. 2007). As has been noted before, comets are thought to contain some 
of the most pristine material from the early solar system, and through the study of comets, 
we can try to determine the history and development of the solar system. 

In the beginning, the expected target for the mission was comet 46P/Wirtanen 
(Schwehm and Schulz 1999). Comet 46P/Wirtanen is thought to be a small but highly 
active comet, with layers of exposed material on its surface as a result of the complete 
removal of material from the strong outgassing (Groussin and Lamy 2003). However, 
due to the failure of an Ariane rocket before take off, the original launch date in 2003 
was delayed until March 2004 (Glassmeier et al. 2007). Consequently, the delay meant 
that 46P/Wirtanen was no longer on a favourable orbit and a new target comet had to be 
found. In this case, comet 67P/Churuyumov-Gerasimenko was chosen, and this comet is 
discussed in the next section. 
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Despite the change in target, the goals of the Rosetta mission remained unchanged 
since the initial proposal, and are described in Schwehm and Schulz (1999), Glassmeier 
et al. (2007) and Schulz (2009). The mission aimed to: 


globally characterise the comet in terms of dynamic properties, surface morphology 
and composition 


determine the chemical, isotopic and mineralogical composition of volatile and re- 
fractory material in the nucleus 


determine the physical properties of the volatile and refractory material, and their 
interrelations 


probe the development of cometary activity and surface processes in the nucleus 
and the inner coma 


study the evolution of outgassing and solar wind interactions during perihelion voy- 
age 


characterise main belt asteroids in terms of dynamic properties, surface morphology 
and composition 


The Rosetta spacecraft was finally launched on March 2, 2004 from Kourou, French 
Guiana, and consisted of an orbiter, Rosetta, and a lander, Philae. The orbiter contained 
eleven instruments (Glassmeier et al. 2007), as shown in Figure 1.4, including: 


Alice (ultraviolet imaging spectrometer) 

CONSERT (COmet Nucleus Sounding Experiment by Radio wave Transmission) 
COSIMA (COmetary Secondary Ion Mass Analyser) 

GIADA (Grain Impact Analyser and Dust Accumulator) 

MIDAS (Micro Imagaing Dust Analysis System) 

MIRO (Microwave Instrument for the Rosetta Orbiter) 

OSIRIS (Optical, Spectroscopic and Infrared Remote Imaging System) 

ROSINA (Rosetta Orbiter Spectrometer for Ion and Neutral Analysis) 

RPC (Rosetta Plasma Consortium) 

RSI (Radio Science Investigation) 


VIRTIS (Visible and InfraRed Thermal Imaging Spectrometer) 


In addition, the Philae lander housed a further ten instruments: 


APXS (Alpha Proton X-ray Spectrometer) 
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Figure 1.4: Schematic diagram of the Rosetta orbiter with instruments labelled (Credit: 
Pline CC BY-SA 4.0). 


e CIVA (Comet Infrared and Visible Analyser) 

e CONSERT (COmet Nucleus Sounding Experiment by Radio wave Trasmission) 
e COSAC (COmetary SAmpling and Composition Experiment) 

e MUPUS (MUlti-PUrpose Sensors for surface and sub-surface science) 

e Ptolemy (geochemistry of light elements) 

e ROLIS (ROsetta Lander Imaging System) 

e ROMAP (ROsetta lander Magnetometer And Plasma monitor) 

e SD2 (Sampling, Drilling and Distribution) 

e SESAME (Surface Electric Sounding and Acoustic Monitoring Experiment) 


Between 2005 and 2010, Rosetta performed several gravity assists around Earth and 
Mars, and flew passed two asteroids: on September 5, 2008, Rosetta flew by aster- 
oid (2867) Steins (Glassmeier et al. 2007) and on July 10, 2010 around (21) Lutetia 
(Thomas et al. 2012). The next stop for Rosetta was the main target of the mission, comet 
67P/Churyumov-Gerasimenko. After being in hibernation for three years, the spacecraft 
woke up in January 2014 and performed rendezvous manoeuvres at 4.5-4.0 AU distances 
in preparation for the start of the near-nucleus phase, which began in August 2014 (Glass- 
meier et al. 2007). The first Rosetta results come from this time in July-August 2014, 
some of which are summarised in Section 1.4 and 1.6. 

On November 12, 2014, Rosetta delivered Philae to the comet surface. However, this 
part of the mission was not entirely successful. Due to a failure of the anchor harpoons, 
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Figure 1.5: Image from the CIVA instrument on Philae of the comet surface (Credit: 
ESA/Rosetta/Philae/CIVA). 


which did not fire, and the cold gas system, which could not push the lander onto the 
surface, Philae did not land as intended onto the surface and bounced several times (Biele 
et al. 2015). For a long time, Philae was thought to be lost until it was finally spotted 
in September 2016, lying on its sight on an uneven part of the comet surface (Schröder 
et al. 2017). Despite the suboptimal landing of Philae, it still managed to transmit many 
interesting close-up images of the nucleus as shown in Figure 1.5. 

From November 2014 until September 2016, Rosetta escorted the comet along its 
perihelion passage, getting as close as 1.24 AU from the sun on August 13, 2015. The 
nominal end of the mission was supposed to be December 2015 (Glassmeier et al. 2007) 
but the performance of Rosetta allowed the mission to continue into 2016. The mission 
ended on September 30, 2016 with a controlled impact onto the nucleus. 

Over the two year mission, Rosetta examined comet 67P/Churyumov-Gerasimenko in 
unprecedented detail and over an extended period of time. Never before had a comet been 
observed in such an intensive manner, and during the mission, some of its mysteries were 
revealed. 


1.4 67P/Churyumov-Gerasimenko 


Comet 67P/Churyumov-Gerasimenko was first discovered in 1969 by Klim Ivanovich 
Churyumov and Svetlana Ivanova Gerasimenko (Lamy et al. 2007). Whilst working at 
the National Taras Shevchenko University of Kiev, Churyumov and Gerasimenko trav- 
elled to the Astorphysical Institute in Almaty, Kazakhstan to make observations of sev- 
eral short period and new comets. On Spetember 11, 1969, Gerasimenko attempted to 
observe comet 32P/Comas Solà but the photographic plate from their observation seemed 
to be defective. Nevertheless, an object in the frame appeared to be 32P/Comas Solà, so 
Gerasimenko kept the plate for analysis back in Kiev. 

When they returned in October and studied their observations, Churyumov and Gerasi- 
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Figure 1.6: NavCam image of comet 67P/Churyumov-Gerasimenko from March 2, 2015 
whilst Rosetta was at a distance of 90 km from the comet centre. The image has been 
enhanced to show the activity of the comet (Photo credit: ESA/Rosetta/NAVCAM, CC 
BY-SA IGO 3.0). 


menko found that there were actually two comets in the plate from September 11: 
32P/Comas Solà and another unidentified comet. They had inadvertently made four fur- 
ther observations of this new comet between September 9-21 and this new discovery was 
reported to the International Astronomical Union. The first sighting of this comet was then 
verified by John Bortle after studying an image from Charles Scovil taken on September 
14, 1969, and the new comet became 67P/Churyumov-Gerasimenko. 


Due to a favourable orbit, comet 67P had been seen as a potentially accessible object 
for a long time (Yeomans and Chodas 1989) but still relatively little was known about it 
at the time. 

The orbit of comet 67P/Churyumov-Gerasimenko was found to have quite a chaotic 
history as a result of repeated encounters with Jupiter, and the last significant planetary 
encounter changed the orbit quite dramatically (Krolikowska 2003, Maquet 2015). The 
perihelion distance decreased from 2.74 to 1.28, the eccentricity increased from 0.36 to 
0.63, and the orbital period decreased from 9 to 6.6 years. The orbit is reasonably well 
defined to 200 years into the future and 200 years into the past, but beyond that, encounters 
with the giant planets make it difficult to reconstruct the exact orbit. 
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Figure 1.7: The shape model of comet 67P/Churyumov-Gerasimenko as found by 
(Preusker et al. 2017). The red line denotes the rotation axis which is tilted at 52? and 
defines the latitude-longitude co-ordinate system on the surface. 


Comet 67P/Churyumov-Gerasimenko attracted much more interest in 2003 when it 
was announced as the target comet for the Rosetta mission. Between the initial observa- 
tions in 1969 and the announcement in 2003, the comet was only observed in three cam- 
paigns in 34 years but observation time was then found in 2003, 2004 and 2005 (Lamy 
et al. 2007). From initial lightcurve data taken by the Hubble Space Telescope in 2003, 
the shape of 67P/Churyumov-Gerasimenko was determined and found to be non-spherical 
(Lamy et al. 2006), the first indication of the unusual shape of comet 67P/Churyumov- 
Gerasimenko. As a result of the mission, much more is known about the comet than any 
other small body in our solar system. 


From the OSIRIS and navigational camera images, the detailed shape of comet 
67P/Churyumov-Gerasimenko can be seen (Figure 1.6). By analysing over 1500 images, 
Preusker et al. (2017) used stereo-photogrammetric to derive a high resolution 3D model 
of the comet. Figure 1.7 shows the result of this work, although only a decimated version 
with 125,000 triangular is shown, rather than the full 1 m scale model with 44 million 
facets. 

One of the most striking aspects of 67P is that it is a bilobed comet with a small 
head lobe and large body lobe which is in stark contrast to the more spherical comets 
previously observed in Figure 1.3. The origins of this shape are still unknown. Some 
argue that the bilobed comet is primordial as a result of low velocity accretion during the 
early formation of the solar system (Massironi et al. 2015, Jutzi and Asphaug 2015), but 
others say that based on its size, 67P would have experienced catastrophic collisions since 
then and that the shape need not be pristine (Morbidelli and Rickman 2015). The bilobed 
shape may then come from collisional disruptions later in the history of the solar system 
(Rickman et al. 2015, Jutzi and Benz 2017). This is still a matter of debate. 
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Figure 1.8: The regions on the surface of comet 67P/Churyumov-Gerasimenko as defined 
by (El-Maarry et al. 2015). These are: (a) Hapi, (b) Seth, (c) Ash, (d) Ma'at, (e) Babi, 
(f) Aten, (g) Hathor, (h) Serqet, (1) Nut, (j) Anuket, (k) Hatmehit, (1) Apis, (m) Atum, (n) 
Bastet, (0) Maftet, (p) Neith, (q) Aker, (r) Khepry, (s) Anubis, (t) Wosret, (u) Imhotep, (v) 
Sobek, (w) Khonsu, (x) Anhur, (y) Geb, (z) Bes. 


OSIRIS also revealed the geological variety of the nucleus. As documented in El- 
Maarry et al. (2015) and El-Maarry et al. (2016), the surface of comet 67P/Churyumov- 
Gerasimenko can be divided into 26 distinct morphological regions based on terrain and 
structure (Figure 1.8). For example, the Imhotep region is one of the most geologically 
diverse, with a smooth bouldered surface as well as rougher terrain, evidence of terracing, 
and filled and unfilled circular features. In contrast, Hathor on the head lobe is cliff-like, 
and has aligned lineaments and fractures. There are dust covered terrains (Ma'at and 
Ash regions) and regions of brittle consolidated material (Seth) (Thomas et al. 2015). 
Numerous features are apparent including depressions, wind tails and dune-like features, 
suggesting that a variety of processes occur on the comet surface. 


The rotation axis of 67P is marked on Figure 1.7 with a red line. The poles are located 
in the neck region between the head and body lobes. From this axis, lines of longitude 
and latitude are defined as for Earth, allowing the map in Figure 1.8 to be made. However, 
due to the complex shape of 67P, the lines of longitude and latitude are quite irregular. 
The equator runs around the body lobe, through the neck region and around the head lobe. 
The rotation axis is tilted at an angle of approximately 52° of its orbital plane, and with 
such a severe tilt, there are pronounced seasonal effects between the north and south poles 
which will be discussed in Chapters 3 and 4. 


The ROSINA spectrometer has also provided a wealth of information regarding the 
composition of the coma. ROSINA has the capabilities to measure the abundance of a 


27 


1 Introduction 





host of oxygen, nitrogen and sulphur bearing molecules (Le Roy et al. 2015), and one of 
the key findings points towards a dichotomy in the CO and CO; activity, with significant 
enhancements in these molecules from the southern hemisphere rather than the northern 
hemisphere. The ROSINA data also enabled the time variability of the abundance for 
a variety of species to be tracked (Hansen et al. 2016, Hässig et al. 2015), and the D/H 
(deuterium/hydrogen) ratio to be determined, being 3 times higher than the terrestrial 
value and higher than that found for other Jupiter Family comets (Altwegg et al. 2015). 

In addition, the VIRTIS instrument was able to identify the presence of non-volatile 
organic compounds on the surface of the comet from the observed low reflectance, spec- 
tral slopes, and the absorption feature between 2.9-3.6 um (Capaccioni et al. 2015); 
while De Sanctis et al. (2015) were able to map the VIRTIS temperatures on the sur- 
face and hence derive a diurnal cycle for water to explain the observed outgassing. An- 
other instrument, RSI, was used to determine the bulk mass and density of the comet - 
(9.982 + 0.003) x 10? kg and (533 + 6) kg m^? - from velocity perturbations at fly- 
by distances and using the comet volume measured by Preusker et al. (2015). RSI also 
determined the global gravity field (Pátzold et al. 2016). In addition, this work implied 
that the interior must be fairly homogeneous, as did data from CONSERT (Kofman et al. 
2015), which used the radar signal from Philae to find that the interior should be relatively 
homogeneous on ~10 m scales. 

The dust component of 67P/Churyumov-Gerasimenko has also been characterised 
by various instruments onboard Rosetta and Philae. There is evidence for a dusty crust 
(Schulz et al. 2015) from the GIADA instrument, which also measured the mass distribu- 
tion of dust grains (in the range 10? - 10? kg) and the dust-gas ratio, finding a value of 
4 + 2 (Rotundi et al. 2015). COSIMA was able to collect and analyse the dust particles 
and determined that there are whole typological classes of dust, with many of them in 
clusters, and similar to micrometeorites and interplanetary dust particles (Langevin et al. 
2016). 

The results briefly mentioned here are only a fraction of the findings that the Rosetta 
mission has revealed about comet 67P/Churyumov-Gerasimenko. Much more research 
has been done, and will be done, with the data from these instruments as well as others, 
including MIRO, the Microwave Instrument for the Rosetta Orbiter. 


1.5 The Microwave Instrument for the Rosetta orbiter 


The Microwave Instrument for the Rosetta Orbiter (MIRO), was one of the eleven in- 
struments onboard the Rosetta spacecraft. This section provides a brief overview of the 
instrument, but the full technical details can be found in Gulkis et al. (2007). 

MIRO consisted of a primary 30 cm reflector telescope and two heterodyne receivers 
operating at centre band frequencies of 190 GHz (1.6 mm) and 562 GHz (0.5 mm). These 
are respectively referred to as the millimetre and submillimetre receivers hereafter. Both 
receivers had broadband continuum channels for the measurement of sub surface bright- 
ness temperatures. The brightness temperature is defined as the required temperature of a 
blackbody which fills the MIRO beam to produce the observed power. MIRO measured 
an antenna temperature, 7'4, first, which was calculated from the received power, P, the 
Boltzmann constant, k, and the bandwidth, dv: 
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Table 1.1: Molecular transitions observed by the MIRO submillimetre spectrometer, re- 
produced from Gulkis et al. (2007) 




















Molecule Frequency (GHz) Transition 
Ammonia 
NH; 572.498 J(1-0) 
Carbon Monoxide 
co 576.268 J(5-4) 
Methanol 
CHOH 553.146 8(1)-7(0) E 
CHOH 568.566 3(-2)-2(-1) E 
CHOH 579.151 12(-1)-11(-1) E 
Water 
H,!°O 556.936 1(1,0)-1(0,1) 
H;"O 552.021 1(1,0)-1(0,1) 
H;'?O 547.676 1(1,0)-1(0,1) 
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Figure 1.9: Schematic representation of the MIRO units. There are four units: Sensors 
Unit, Sensor Backend Electronics Unit, Electronics Unit and an Ultra-Stable Oscillator 
(USO) Unit. The Sensors Unit contained the telescope, baseplate and the optical bench 
(the calibration targets, the 3-position mirror, the quasi-optics, and the two heterodyne 
receivers). The Sensor Backend Electronics Unit held the Intermediate Frequency Pro- 
cessor (IFP), phase lock loop, sensor electronics and frequency synthesisers. The Elec- 
tronics Unit contained the Chirp Transform Spectrometer, the instrument computer, power 
circuits and spacecraft interface circuits. The USO Unit was a self contained unit which 
held the oscillators used in the downconversion of the received signal. Image inspired by 


(Gulkis et al. 2007). 
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Figure 1.10: Technical diagram of the heterodyne receivers of MIRO. The front ends of 
the millimetre and the submillimetre receivers both have the same design. The power is 
detected at RF radio frequencies and sent from the optical bench to a subharmonically 
pumped mixer via a feed horn. At the mixer, the RF signal is multiplied by a signal (LO) 
from a low frequency oscillator. The LO frequency is produced by a Gunn oscillator. 
The mixer produces an IF Intermediate Frequency which is passed through a low noise 
amplifier to an IF processor. All of this is controlled by the power unit. Image inspired by 
(Gulkis et al. 2007). 


P 


Tı = 
^ kdv 


(1.1) 


The antenna temperature can then be converted to a brightness temperature, Tg, using: 





1 h 
Tao) = -—— — (1.2) 
Kin * 1) 


where v is the frequency and Ah is the Planck constant. 

In addition, a Chirp Transform Spectrometer (hereafter, CTS, Hartogh and Hartmann 
(1990)) was attached to the submillimetre receiver for the detection of molecular spectral 
lines. The CTS was capable of detecting absorption and emission lines from some of the 
most abundant species previously observed in comets, including: carbon monoxide (CO); 
methanol (CH3OH); ammonia (NH3); and three isotopologues of water (H,!°O, H,!’O, 
H,!8O). The transitions for these molecular lines are listed in Table 1.1. 

MIRO consisted of four units (Figure 1.9): Sensors Unit, Sensor Backend Electronics 
Unit, Electronics Unit and an Ultra-Stable Oscillator Unit (Gulkis et al. 2007). The Sen- 
sors Unit contained the telescope, baseplate and the optical bench. The two heterodyne 
receivers were part of the optical bench, along with the hot and cold calibration targets. 
The Sensor Backend Electronics Unit held the Intermediate Frequency Processor (IFP), 
phase lock loop, housekeeping circuits and frequency synthesisers. The CTS was part 


30 


1.5 The Microwave Instrument for the Rosetta orbiter 





of the Electronics Unit, as well as the instrument computer, power circuits and spacecraft 
interface circuits. Finally, the Ultra-Stable Oscillator Unit was a self contained unit which 
held the oscillators used in the downconversion of the received signal. Figure 1.9 shows 
how the signal flows from the telescope, via the heterodyne receivers, to the IFP and the 
spectrometer. 

The primary 30 cm offset parabolic reflector telescope had a focal length of 37.5 cm 
and formed the external part of the Sensors Unit, with the optical bench mounted inside 
the spacecraft. The telescope was located on the top of Rosetta and aligned with the pay- 
load line of sight. The angular resolution for the millimetre and submillimetre receivers 
were 23.8 arcmin (6.9 mrad) and 7.5 arcmin (2.2 mrad), respectively. At an altitude of 
30 km from the nucleus, the beams would have sizes of 200 m in the millimetre channel 
and 70 m in the submillimetre channel. Light from the primary telescope was reflected in- 
side to the optical bench. At the optical bench, the incoming light was split into millimetre 
and submillimetre components and directed to two heterodyne receivers. 

Heterodyne receivers work by mixing (multiplying) the received frequency, called RF 
(radio frequency), with the frequency of a local oscillator (LO). The multiplication of 
the RF and LO frequencies provides sum and difference frequencies called intermediate 
frequencies (IF). Usually, the IF is used for further processing and it can be chosen in a 
way that it appears in a range where electronic processing is feasible (Janssen 1994). In 
the case of MIRO, the RF needed to be downconverted to lower frequencies. There are a 
number of reasons why it can be beneficial to shift the frequency range lower. Circuitry 
performance is relatively poor at high frequencies and transistors provide little amplifica- 
tion. 

Both the millimetre and submillimetre heterodyne receivers worked in the same way 
and are shown schematically in Figure 1.10. The RF signal arrived from the optics unit 
and went into a subharmonically pumped mixer via a feed horn. At the mixer, the RF 
signal was multiplied by the LO signal to downconvert the signal to the IF. The LO was 
provided by a Gunn oscillator which operated in a frequency switching mode. In this 
mode, the Gunn oscillator frequency switched by +5 MHz every 5 seconds, creating shifts 
in the spectral lines. The shifts could be used to remove gain and baseline variations which 
occur in the receiver and thus improve the sensitivity of the instrument. The IF signal was 
then amplified and finally sent to the IFP for detection. The millimetre and submillimetre 
RFs were in the range 186.7-189.7 GHz and 547.6-579.2 Hz, respectively, whilst the 
LOs were at ~95 GHz and ~282GHz (Gulkis et al. 2007). The IFs then had ranges of 
1-1.5 GHz and 5.5-16.5 GHz, respectively. 

The IFP had two continuum channels, one for each receiver, and prepared the spectral 
bands for the CTS (Hartogh and Hartmann 1990). The CTS was a 4096 channel spec- 
trometer with a 180 MHz bandwidth and a spectral resolution of 44 kHz (v/Av ~ 10’). 
The spectrometer enabled the detection of the eight molecular transitions listed in Table 
1.1. The lines are separated by a filterbank at the input of the CTS. Table 1.2 details 
how the RF was downconverted to the bandwidth of the CTS through a series of LO. 
The downconversion to IF1 is performed in the submillimetre receiver. The remaining 
IFs are produced in the IFP, which uses nine mixers and three LO sources to complete 
the downconversion. The bandwidth of the RF reduces from 31.6 GHz to 180 MHz in 
the CTS. The MIRO CTS spectrometer was uniquely capable of making high resolution 
observations of the intrinsic narrow Doppler broadened spectral lines which result from 
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Table 1.2: Downconversion of the radio frequency RF (by multiplication with local oscil- 
lator frequencies LO) for each molecular transition given in Table 1.1. The RF is down- 
converted to a series of intermediate frequencies IF and finally to the bandwidth of the 
CTS (Hartogh and Hartmann 1990). All frequencies given in units of MHz. Reproduced 
from Gulkis et al. (2007) 








H>!°O | H»'°O | CH;OH | H5*O CO | CHOH] NH; | CHOH 
RF 556936 | 552022 | 579151 | 547676 | 576268 | 568566 | 572498 | 553146 
IF 1 (LO 1) 5877 | 10791 | 16338 | 15137 | 13455 | 5753 9685 9667 
IF 2 (2xLO2) = 6427 : € 9091 1389 | 14049 | 14031 
IF 3 (LO4) z = - - 1363 - 6321 6303 
IF 4 (LO4) x s : : à Z 1407 1425 
IF 5 (LO3) - 1301 i 3 5 - - s 
IF 6 (LO2) 1270 z 9191 7990 4 - - - 
IF 7 (LO3) z - - 5808 - - - - 
IF 8 2xLO2) - = 2044 1339 à - - - 
IF 9 (LO4) - 4 6408 = - 4 4 ^ 
IF 10 : : 1320 2 = - - g 
output to CTS | 1270 1301 1320 1339 1363 1389 1407 1425 



































the low gas temperatures and low pressures around comets. 

The large dataset provided by MIRO allows measurements of molecular column den- 
sities and production rates to be made from the observed spectral lines, as well as the ki- 
netic properties within the inner coma. The MIRO brightness temperatures derived from 
the continuum channel also enable the thermal properties of the subsurface to be inferred. 
The temporal and spatial capabilities of the instrument allow these measurements to be 
taken across the whole surface of the nucleus and throughout the entire two year mission. 


1.6 MIRO aims, results and spectra 


Before the beginning of the Rosetta mission, several measurement objectives were out- 
lined for the MIRO instrument in Gulkis et al. (2007). One objective was to measure the 
abundances of some of the major volatile species detectable by MIRO. The sublimation of 
H20 and CO are thought to be two of the main drivers of activity in comets and during the 
lifetime of Rosetta, the change in abundance of these two key volatiles can be measured 
as 67P/Churyumov-Gerasimenko approaches and recedes from the sun after perihelion. 

Another important objective related to the nucleus itself, as the MIRO data can help 
to improve our understanding of the processes which control cometary activity in the sub- 
surface layer, and hence characterise the surface over the first few centimetres. Knowl- 
edge of the thermal gradient in the upper layers of the nucleus can help to constrain the 
physical nature of the comet and determine how much material could still be pristine from 
the early solar system. 

Furthermore, MIRO offered a unique opportunity to study the evolution and develop- 
ment of the inner coma. This region is impossible to probe from ground-based telescopes 
or in situ instruments, but as a remote sensing unit, the MIRO measurements can be used 
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to estimate the kinetic temperature, expansion velocity and molecular density along its 
line of sight. This will enable the structure and dynamics of the inner coma to be studied 
in great detail. 

The final main objective was to take advantage of the fly-bys of asteroids Lutetia and 
Steins to probe for low levels of gas in these objects. This was accomplished by Gulkis 
et al. (2010) and Gulkis et al. (2012). 

Much of the early MIRO work regarding 67P went into characterising the surface by 
looking at the subsurface temperatures. The submillimetre and millimetre emissions are 
estimated to come from depths of 1 cm and 4 cm respectively (Schloerb et al. 2015). 
Alternatively, Choukroun et al. (2015) propose that the depths of these brightness temper- 
atures may originate from much deeper layers if ice is present in the surface, which could 
extend the electrical penetration depth, giving rise to potential depths of 15 - 20 cm in the 
sub-millimetre channel and 20 - 30 cm in the millimetre channel for icy regions. Schloerb 
et al. (2015) and Gulkis et al. (2015) were able to demonstrate the diurnal variation of the 
subsurface temperatures when 67P/Churyumov-Gerarsimenko was more than 3 AU away 
from the sun. The variation appears to be a function of latitude, owing to the fact that 
the received solar flux strongly depends on the latitude, and there is a possible tempera- 
ture lag, with the highest subsurface temperatures occurring in the early or mid afternoon, 
local time. Night time subsurface temperature measurements are given in Choukroun 
et al. (2015), showing the surface to be very cold, between 17-50 K across both instru- 
ment channels. The low sub-surface temperatures can only be measured precisely with a 
millimetre/sub-millimetre instrument like MIRO operating at this wavelength range. 

Schloerb et al. (2015), Choukroun et al. (2015) and Gulkis et al. (2015) made esti- 
mates for the thermal inertia of the surface, the ability of the surface to resist changes 
in temperature. Materials with a high thermal inertia, such as rocks, change temperature 
slowly, whereas materials with a low thermal inertia, like sand, thermalise quickly. A 
thermal inertia in the range 10 - 30 JK ! m?s 9? for the surface is found by Schloerb 
et al. (2015), 10 - 50 JK^! m?s 9? from Gulkis et al. (2015), and from Choukroun et al. 
(2015), 10 - 40 JK^! m?s^9? in the sub-millimetre channel and 20 - 60 JK^! m^?s-?? in the 
millimetre channel. The discrepancy between the two channels is found when Choukroun 
et al. (2015) try to fit their model to the millimetre and sub-millimetre measurements in- 
dependently, thus obtaining slightly different thermal inertia ranges. They suggest that 
the presence of ice in the sub-surface, which has a wavelength dependent electrical pene- 
tration depth, could explain the difference. 

As well as the nucleus, the MIRO data also offers an opportunity to probe the coma 
around the comet. Using parametrised coma profiles of temperature, velocity and density, 
Lee et al. (2015) measured the outgassing rate, the terminal gas expansion velocity and 
the terminal gas temperature across the nucleus from measurements made in August 2014. 
Interestingly, they found that the outgassing of the comet was not strongly correlated to 
the solar illumination or the surface temperature, implying that illumination is not the 
only driver of cometary activity. It is revealing that Lee et al. (2015) report an increase 
in the outgassing just after sunset on the comet, suggesting that the deposits of ice must 
be below the surface but within the penetration depth of the thermal wave. In addition, 
they find that the neck is the most active region at this time, similar to findings by other 
instruments. The spatial variation results in some areas which have 30 times more water 
production than less active areas. The temporal variation changes by only a factor of 
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Figure 1.11: Spectral lines observed by MIRO of (a) H,!°O, (b) H,'’O, and (c) H,'?O, 
on July 1, 2015 at 04:45 with a 30 minute integration time. Doppler shift velocity is 
depicted on the x axis and antenna temperature on the y axis. The measured continuum 
temperature for these observations was 173 K. The viewing geometry is shown in (d) with 
the red line and circle indicating the MIRO beam and the MIRO footprint, respectively. 
The yellow line is the vector towards the sun. 


2 though, so the spatial distribution is much more prevalent than the temporal changes 
during this early study period. 


Biver et al. (2015) also derived the column densities around the nucleus in the early 
part of the mission from the MIRO spectra. They are able to map the the spectra to 
the pointing of the spacecraft and thus create a map of temperature, expansion velocity 
and water column densities around comet 67P/Churyumov-Gerasimenko. The results are 
similar to Lee et al. (2015), in that the neck region is found to be the origins of higher 
column densities than other parts of the comet. There is also evidence for low levels of 
production on the night side, although this could also be a result of back flow towards the 
dark unilluminated side. 


Figure 1.11 shows examples of the spectral lines observed by MIRO of three water 
isotopologues from observations at 04:45 on July 1, 2015 with an integration time of 30 
minutes. H5'*O and H5"O are minor species of water with terrestrial abundance ratios 


of and —1., respectively, relative to H,!°O (Baertschi 1976, Jabeen and Kusakabe 


I 
498.7 2632? 
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Figure 1.12: Spectral lines and viewing geometry of observations from MIRO on July 
1, 2015 with 30 minute integration time. (a) H5!^O in absorption against warm nucleus 
from 04:30-05:00, subsurface temperature of 173 K and viewing geometry shown in (b) 
which is the same as in Figure 1.11. (c) H2!°O in emission against cold space from 14:00- 
14:30, background temperature ~5 K and viewing geometry shown in (d). (e) H2!*O in 
absorption and emission MIRO observes nadir and limb simultaneously as the viewing 
geometry shows in (f). The viewing geometry shows a red line and circle indicating the 
MIRO beam and the MIRO footprint, respectively. The yellow line is the vector towards 
the sun. 
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1997). They are shown here in absorption against the warm nucleus surface. All of 
these molecules consist of two hydrogen atoms and an oxygen atom but they differ in 
the number of neutrons within the oxygen atom, denoted by the superscript number. The 
H,!°O line, the main species of water is shown in Figure 1.11a, and the wide deep shape 
of the spectral line shows that it is saturated and optically thick. The H'*O line shown 
in Figure 1.11c is optically thinner, thus narrower and having a lower contrast than the 
H,!°O line. The spectrum in Figure 1.11b shows the hyperfine structure in the H,'7O 
line. It is five times less abundant than the H,!*O and the line is weaker by this amount 
(requiring 25 times the observing time). Figure 1.11d shows the viewing geometry of this 
observation, which at this time was in a terminator orbit, observing close to the transition 
from day to night side. 

The MIRO spectral lines can be observed in absorption and emission, as seen in Figure 
1.12 for the optically thick H5!^O line. In Figure 1.12a, MIRO makes a nadir observation 
and the line is seen in absorption against the warm nucleus with a temperature of 173 K on 
July 1, 2015 at 04:30. The viewing geometry is shown in Figure 1.12b and is for the same 
time as in Figure 1.11. Due to the high optical depth (> 10), the line core is saturated, 
and the velocity profiles makes the line shape asymmetric. There is also evidence of an 
emission wing at ~1 km/s. 

The H5/5O line in emission from a limb observation is shown in Figure 1.12c with 
the viewing geometry in Figure 1.12d. The emission line is quite complex owing to the 
optical thickness and the velocity profile. The high opacity line creates the double peaked 
features, where an optically thin line would be approximately Gaussian in shape. The 
non-constant velocity profile then skews the emission line, making it asymmetric. The 
line core in this example occurs at ~ -1 km/s. 

Figure 1.12f shows MIRO observing partly on and partly off the nucleus, and as a re- 
sult, the spectra in this time period exhibits both absorption and emission features (Figure 
1.12e) since it is a convolution of nadir and limb measurements. 


1.7 Thesis aims 


In this work, I explore aspects of the first three aims outlined in the previous section. In 
Chapter 3, the abundance of water is measured by tracking the change in the line area 
ratio of the H,'°O and H3!?O lines to find spatially resolved local water production rates. 
This is done from ~3 AU pre-perihelion to ~2 AU post-perihelion, allowing the change in 
production rate to be seen during the comets close approach to the sun. The meaning of 
water production rate curves and their derived slopes is investigated in Chapter 4. Chapter 
5 describes how the thermal gradient can be found from the MIRO brightness temperature 
measurements, and hence derive estimates for the thermal inertia of the comet surface. In 
the chapter on Future Work, the derivation of temperature, velocity and density profiles in 
the inner coma is discussed, showing how the spectra can be inverted to determine these 
parameters. Finally, in Chapter 7, the results presented in the previous chapters are briefly 
summarised in their wider context. 

The next chapter contains a brief introduction to a few scientific concepts which will 
occur in later chapters. 
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In this chapter, a few basic physical and theoretical concepts will be introduced which 
will be mentioned in future chapters. 


2.1 Haser model 


The Haser model is the simplest and most widely used model to describe the outflow 
of gas from a comet nucleus. It assumes a spherical comet has a global production rate 
for parent molecular species and that there is an isotropic outflow of material with finite 
lifetimes for the gas. The density distribution of gas produced from the nucleus is given 
by the equation: 





N(r) = exp(-r/R) (2.1) 


Anr?v 

where N(r) is the density at a distance, r, Q is the production rate of the gas in 
[molecules s~'], v is the expansion velocity of the gas into the coma, and R is the scale 
length of the gas, given by R = vr where r is the lifetime of the molecule (De Pater 
and Lissauer 2015). Here, the velocity profile is assumed to be constant and move ra- 
dially outwards, although the velocity can also depend on the distance from the comet. 
For molecules with a long lifetime and thus long scale length, the density yields a 1/7? 
dependence when following the assumptions given above. 

The density distribution of daughter species can also be calculated with the Haser 
model if the lifetimes and production pathways of these daughter species are known as 
well. 

The Haser model is thought to be a reasonable first approximation when deriving the 
production rate from comets (Snodgrass et al. 2017), and it has been reasonably success- 
ful in the past for deriving the density around a nucleus (Cochran 1985). However, the 
assumptions stipulated above, such as sphericity and radial outflow, can lead to problem- 
atic attempts to fit observed spectra with the Haser model formalism and other methods 
can be more suitable (Bockelée-Morvan and Gérard 1984). 


2.4 Atwo level atom 
An atom with two energy levels is represented in Figure 2.1. It has a lower level / with an 
energy E, and an upper level u with an energy E,,. Photons can be absorbed and emitted 


via transitions between the two energy levels. There are electronic transitions, vibrational 
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Figure 2.1: Schematic representation of an energy level system in an atom. Stimulated 
absorption - a background radiation field inducing a transition from the lower, /, level 
to the upper, u, level through the absorption of a photon with energy Avo 2 E, — Ej. 
Stimulated emission - a background radiation field inducing a transition from the upper 
level to the lower level resulting in the emission of a photon with energy Avo — E, — Ej. 
Spontaneous emission - a spontaneous transition takes place between the upper and lower 
energy levels resulting in the emission of a photon. Collisional excitation - collisions 
between atoms result in a transition from the lower state to the upper state. Collisional 
de-excitation - collisions between atoms result in a transition from the upper state to the 
lower state. 
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transitions and rotational transitions. Rotational transitions are particularly important for 
the microwave regime. 

The first case demonstrates stimulated absorption, whereby the energy from an in- 
coming photon is absorbed by the atom and there is a transition from the lower level to 
the upper level. This is controlled by the Einstein coefficient B;, [J s ! m? sr molec'!] 
which is the probability per second per spectral energy density of the radiation field that 
a transition will occur. For stimulated absorption, the absorbed photon has an energy of 
hvo 2 E,— Ej, where vo is the frequency of the transition [Hz] and Ah is the Planck constant 
[6.626x10^7^ J Hz! ]. 

Stimulated emission is controlled by the Einstein coefficient B,; [J s ! m? sr molec !], 
giving the probability per second per spectral energy density of the radiation field that a 
transition will take place from the upper level to the lower level. 

Spontaneous emission is described by the Einstein coefficient A,; [s ! molec !] which 
is the probability per second that a transition will spontaneously occur from the upper 
level to the lower level. For stimulated and spontaneous emission, the emitted photon has 
an energy of Avo — E, — E, and the inducing radiation field for stimulated emission must 
have a frequency close to vo. 

In addition to these three radiative transitions, transitions between the lines can also 
be induced by collisions between atoms. In this case, the transitions between the lower 
and upper energy levels, and the upper and lower levels are represented by collisional 
rates C and Cu, respectively, which describe the probability per second of a transition 
occurring as a result of a collision. 

The population of individual levels must be solved using a system of equations (Ya- 
mada et al. 2018): 


dn; - = 
a 0= 2 uA; - (niBi; - njBj)J ji] — 2 más - (njB i - niBij)Jij] * 26i -njCi;) 
pi J<i J 
(2.2) 
where n; is the population of energy level i, t is time [s], and J;; is the integrated mean 
intensity: 
1 co Art 
Ji; = — d 1,0;;,d€2 2.3 
S - 


where ¢j;;, is the frequency dependent line emission profile. Equation 2.2 is called the 
statistical equilibrium equation (Rybicki and Hummer 1991). 


2.3 Optical depth 


When an incoming ray of light moves through a medium (Blundell and Blundell 2009), 
the intensity of the electromagnetic wave is changed from its initial spectral intensity, 7,0, 
to its outgoing spectral intensity, /,, in units of [J s ! m? sr ! Hz !]. The subscript , de- 
notes the frequency dependence of the terms. The amount of change in intensity depends 
upon the mass density of absorbing molecules in the material, p, [kg/m?]; and the com- 
position of the medium. The composition of the medium is represented by a parameter 
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called the opacity, x, [m?/kg], which describes how efficiently radiation is absorbed by 
the material. It is also frequency dependent. 

In a homogeneous medium with absorption only, the amount of intensity lost in the 
medium is represented by: 


dl, = -k,pal,ds (2.4) 


where ds is some spatial element of the medium. Integrating Equation 2.4 gives: 


I, = v,0 exp(—k,paS) (2.5) 


which describes an exponential decay of the intensity with distance. This result is 
commonly known as the Beer-Lambert law. The mean free path, /, can be defined as the 
distance for the intensity to decrease by 1/e or: 


1 


KyPa 


l= 





(2.6) 


in the homogeneous medium. 
Optical depth, t,, combines opacity and density into one parameter and defines the 
ability of the material to block light: 


Be «pads | a,(s) ds (2.7) 
0 0 


where s, is the length of the medium and a,(s) is the absorption coefficient [m~']. For 
a non-homogeneous medium, k,, 0, and a, change within the material. 

Notice that the optical depth is a function of the opacity, so that as the opacity in- 
creases, so does the optical depth. However, a material can have a large opacity but a 
small optical depth, depending on the magnitude of the density. 

There are a few special cases for r,. When r, « 1, the medium is said to be optically 
thin and essentially transparent as the medium does not absorb the incoming intensity. In 
this case, the radiation mean free path is far greater than the length of the medium and 
radiation can propagate without any absorption. Any changes in extinction are propor- 
tional to the amount of material in the medium, so in the optically thin regime, extinction 
is proportional to the density of the absorbers. 

When r, = 1, the mean free path is the same length as the medium. Generally, we can 
see about one mean free path, or one optical depth, into a material. 

In the other extreme, where T, » 1, the medium is completely opaque. The medium 
is far greater in length than the mean free path and only radiation from the surface can 
be seen. The proportionality to density is lost in the optically thick case, and the medium 
now acts like a blackbody, radiating with the intensity found at the surface. 

A blackbody is an ideal object, absorbing and emitting with 100% efficiency. Emissiv- 
ity is the ratio of actual emission to blackbody emission. Blackbodies have an emissivity 
equal to 1 and non-blackbodies have an emissivity less than 1. A similar definition exists 
for absorptivity, and Kirchoff’s law states that the emissivity of a medium is equal to its 
absorptivity at a particular frequency. 
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2.4 Radiative transfer equation 


The optical depth is the fundamental property regarding radiative transfer. Radiative trans- 
fer describes the propagation of radiation through a medium (Chance and Martin 2017). 
The change in intensity of radiation as it moves through a medium is modulated by absorp- 
tion, which reduces the intensity, emission, which enhances the intensity, and scattering 
processes, which can reduce and enhance the intensity. The radiative transfer equation 
takes the form: 


2 2 ifs) — asi) (2.8) 
s 

where j,(s) is the spectral emissivity with units of [J s! m? sr ! Hz !] Janssen 1994). 
All of these parameters are functions of frequency, v, and the distance that the radiation 
has moved through the medium, s. The equation describes how the intensity which enters 
the medium is changed by the medium (left hand side) through gains as a result of emis- 
sion from local contributions (j,(s)) and losses from self-absorption (a,(s)/,(s)). In the 
case where emission is negligible, j,(s) = 0, the Beer-Lambert Law in Equation 2.5 can 
be retrieved from the radiative transfer equation. 

The spectral emissivity and the absorption coefficient depend upon the occupation of 
molecular energy levels: 





js) = Fh Au vo Cv, vo) (2.9) 
TT 


"OR eh Bu — f, Bu) hvo Wr, vo) (2.10) 


where f, is the fraction of molecules in the upper energy level and f, is the fraction 
of molecules in the lower energy level, p is the molecular number density [m ?], v is the 
frequency and y is a normalised line profile so that b Wy, vo) dv = 1. The normalised 
line profiles for j,(s) and a,(s) are assumed to be the same (complete frequency redistri- 
bution). The spectral emissivity and the absorption coefficient are functions of the three 
Einstein coefficients from Section 2.2 and they can be expressed in terms of g, and gi, the 
statistical weights of each energy level, using the equations: 


B u u 
Du _ Eu (2.11) 
Bu  & 
and: 
2hy? 
Aul == : By (2.12) 
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where c is the speed of light in a vacuum. 
A source function, $,(s), can be defined as S,(s) 2 so that Equation 2.8 can be 
rewritten as: 


d/,(s) 
dr, 
and the differential is now in terms of dr since dr, = a,(s)ds following from Equation 
2:1. 





- $,(5) - L(s) (2.13) 
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Using the derivation in Appendix A, the radiative transfer in Equation 2.8 can be 
expressed as: 


I, = I,o exp(-T,(8o)) + li 7 Sy exp(-T,(s)) ay(s) ds (2.14) 
0 


The intensity at the observer, /,, is the sum of the attenuated intensity from the source, 
I,o, and all the contributions from the material along the line of sight with a length s; 
(second part of the right hand side of the equation). 

With Equations 2.9 and 2.10, we can express the source function as: 


2h? 1 
Ber c? &fi _ 
$1 Fu 
The population of the energy levels is therefore fundamentally important in order to 
solve the radiative transfer equation. 





(2.15) 


2.5 LTE vs. non-LTE 


A system is in Local Thermodynamic Equilibrium (LTE, hereafter) when the energy level 
distribution is determined by the collisional terms (C;, and C,,) in the statistical equilib- 
rium equation (Equation 2.2). A local region can then equilibrate via collisions and the 
region can be described by a single temperature (Chance and Martin 2017). This may 
occur at a specific altitude in an atmosphere or coma. In LTE conditions, temperature can 
still vary globally within the system but at a specific location, thermodynamic equilibrium 
can be assumed (Blundell and Blundell 2009). 

The assumption of LTE allows the Boltzmann distribution to be applied to the system 
to find the fraction occupying an energy level, i, at a temperature, T: 


-E; 
rar) Ze) (2.16) 
: ZT) 
where Z(T) is the partition function: 
ES 
ZT) = 25 gi exp( =") (2.17) 


which describes how populations are distributed among the different levels and ener- 
gies (Chance and Martin 2017). 
In LTE, the source function can be expressed as: 


2h? — 1] 
S, = — (2.18) 
C exp( #) -1 


which is simply a Planck function, B,(T), where k is the Boltzmann constant and T is 
the blackbody temperature. 

At low microwave frequencies, and thus long wavelengths (4), the Rayleigh-Jeans 
approximation can be used (Janssen 1994), whereby hv « kT, so that the Planck equation 
reduces to: 
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2y?kT 
B(T) - 5 (2.19) 
c 
or: 
2kcT 
BıT) = 7 (2.20) 
because in the Rayleigh-Jeans limit: 
: SL (2.21) 


exp&)-1 mv 


Brightness temperature, T, is defined as the blackbody temperature required in order 
to produce a measured power and can be expressed as: 


2 


C 
TY) = aah, (2.22) 


which allows the radiative transfer equation to be written in terms of the brightness 
temperature: 


TV) = Th ov) exp(—T(So)) + J : T(s) exp(—t,(s)) a@,(s) ds (2.23) 
0 


This form of the equation is commonly used in LTE microwave remote sensing and in 
applications for analysing the MIRO dataset. 

At low densities however, especially those expected for cometary comae, the colli- 
sional frequency of molecules is not high enough to keep the LTE conditions and other 
processes have to be taken into account. In this case we speak of the breakdown of LTE 
and must instead use the non-Local Thermodynamic Equilibrium (non-LTE, hereafter) 
regime. 

When collisions no longer dominate, Equation 2.16 can not be used in non-LTE and 
instead, we must account for all the detailed processes populating and depopulating the 
levels which are expressed in the statistical equilibrium equation (Rybicki and Hummer 
1991). Clearly, having to use Equation 2.2 for each level makes non-LTE calculations far 
more computationally intensive and time consuming than simply using the LTE approx- 
imation. In LTE, the collisional terms dominated so that Equation 2.2 greatly reduced in 
complexity. 

Within cometary comae, it is vitally important to consider the radiative transfer prob- 
lem in non-LTE due to the tenuous nature of the atmospheres (Yamada et al. 2018). 
Nevertheless, for high production comets (>1x10?’ molec/s), the collisional terms can 
sometimes dominate up to 50-100 km, in which case using LTE conditions is justified. 
In addition, the choice between using LTE or non-LTE can depend on the science ques- 
tion: if considering just the line area, LTE may be valid within some uncertainties (as 
most of the line is formed inside the LTE region), but if the information provided by the 
line shape is investigated, then non-LTE calculations have to be considered (such that the 
region around the line core is accurately considered). 

Scattering processes are also important in the radiative transfer equation but for the 
microwave regime which we are interested in, scattering can be ignored. Figure 2.2 is 
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Figure 2.2: Schematic representation depicting where different scattering regimes domi- 
nate at different wavelengths and for different particle radii. The red lines approximately 
indicate the wavelength of MIRO and the size of dust particles. Adapted from Chance 
and Martin (2017) 


reproduced from Chance and Martin (2017) and shows where different scattering regimes 
dominate for various wavelengths and particle sizes. The MIRO instrument takes mea- 
surements at millimetre and submillimetre wavelengths where scattering is mostly negli- 
gible for the small dust particles expected in cometary comae. 


2.6 Inversion methods 


A forward model can be thought of as using a set of parameters, x, in a function, F(x), to 
predict a set of measurements, y. This can be written as: 


y = F(x) (2.24) 


An inverse problem, as the name suggests, is the reverse procedure of the forward 
model: using the measurements to predict the initial parameters. The inverse problem 
takes the form: 


x-2Fy (2.25) 


The function F is often non-linear and hence requires a linear approximation to sim- 
plify the problem. This approximation is called quadrature (Twomey 2013) and F is recast 
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as the weighting function K. Equations 2.24 and 2.25 hence become: 


y=Kx (2.26) 
x=K'y (2.27) 


An inverse problem is said to be well-posed if a solution exists, which is unique, and 
depends continuously on the data. In reality, inverse problems are rarely well-posed. 

In the inverse problem which will occur later in this work, the parameters which we 
want to find, our x, will be temperature, velocity and density values (discretised in the 
radial direction), while the measurements, our y, will be the intensity of water spectral 
lines. As will be discussed later, the desired parameters cannot all be uniquely determined. 
The source of the non-uniqueness arises due to the nature of the radiative transfer equation 
which relates the desired parameters, and due to measurement uncertainty. 

In addition, the measurements may not depend continuously on the parameters. This 
problem is especially prevalent in the optically thick case where changes in the density do 
not always have an effect on the intensity of the spectra. As a result of these issues, our 
problem is said to be ill-posed (Twomey 2013). 

The non-linearity of the problem must be carefully considered because it may, under 
certain conditions, also lead to instability in the inversion. Simply put, this could occur 
when changes in the parameter space do not produce changes in measurement space. 

So, how can these problems be overcome? 

We can rewrite the inverse problem in terms of the least squared solution, where the 
aim is to minimise the equation: 


min |IKx — y|P (2.28) 


In constrained least squares, an additional, if arbitrary, term is added to the cost func- 
tion and Equation 2.28 becomes: 


min |Kx — yIl* * IIyHIF (2.29) 


where y can take any finite value and H is a smoothing matrix, typically the iden- 
tity matrix or a near-diagonal matrix. Adding this additional constraint is called regular- 
ization. This additional term effectively limits the solutions to a subset of all solutions, 
thereby reducing the degeneracy of the problem and also helping to stabilise the inversion. 

Regularization can be performed in a different way that also incorporates correlation 
among the parameters in the a priori information. In can take the form: 


min [Kx — y)$ (Kx — y) + (x = x4)$ 4G — xj) (2.30) 


where S, and S, are measurement and a priori covariance matrices, and x, is some a 
priori information on the parameter space. The covariance matrices can add smoothness 
to the solutions, and thus remove any wild and discontinuous solutions. Furthermore, 
solutions can now be limited to stay close to the a priori. This type of regularisation 
is typically used in optimal estimation methods of inversion (Rodgers 2000). Optimal 
Estimation is an iterative process which takes the form of: 


Xii = Xa € SK (KSSK! 5,) - FÜ) + K(x — xà) (2.31) 
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where i is the iteration number and T denotes a transposed matrix. 

It is challenging to apply Equation 2.31 to radiative transfer problems for comets. 
Firstly, there is little a priori knowledge regarding cometary comae, especially the inner 
regions within 100 km from the nucleus. Creating a suitable x, is therefore difficult. In 
addition, the covariance in the coma is unknown so S, acts like a free parameter which 
would otherwise be fixed if this was an Earth atmosphere inverse problem. In Chapter 6, I 
will describe out the Optimal Estimation inverse method can be applied to measurements 
made from 67P and how these issues can be dealt with. 
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comet 
67P/Churyumov-Gerasimenko from 
the MIRO instrument on Rosetta 


This chapter has been published in A&A as the journal article: Marshall, D. W., et al. 
"Spatially resolved evolution of the local H2O production rates of comet 67P/Churyumov- 
Gerasimenko from the MIRO instrument on Rosetta." Astronomy & Astrophysics 603 
A87 2017. Reproduced with permsission © ESO 


3.1 Summary 


Using spectroscopic and continuum data measured by the MIRO instrument on board 
Rosetta of comet 67P/Churyumov-Gerasimenko, it is possible to derive and track the 
change in the water production rate, to learn how the outgassing evolves with heliocentric 
distance. The MIRO data are well suited to investigate the evolution of 67P, in unprece- 
dented spatial and temporal detail. To obtain estimates of the local effective Haser pro- 
duction rates we developed an efficient and reliable retrieval approach with precalculated 
lookup tables. We employed line area ratios (H2!°O/H»'8O) from pure nadir observations 
as the key variable, along with the Doppler shift velocity, and continuum temperature. 
This method was applied to the MIRO data from August 2014 until April 2016. Perihe- 
lion occurred on August 13, 2015 when the comet was 1.24 AU from the Sun. From the 
start of the period until seven days before perihelion the water production rates increased 
by an order of magnitude, and the derived maximum (1.42+0.51)x10*8 is for August 29, 
2015. The data also indicate that there is an offset in the peak outgassing, occurring 22- 
46 days after perihelion, depending which kind of fitting is used. In the pre-perihelion, 
the production rate changes with heliocentric distance like r,***°, while post-perihelion 
the dependence is r;^?*??. In addition, this method provides well sampled dataset to de- 
termine the spatial distribution of the outgassing versus heliocentric distance. The time 
evolution is definitely not uniform across the surface. During the approach toward the 
perihelion the surface temperature on the southern hemisphere increases rapidly and so 
does the sublimation rate with an exponent of ~ -6 pre- and post-perihelion. We also 
present more detailed regional variation in the outgassing rate, which demonstrates that 
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the highest derived production rates originate from the Wosret, Neith and Bes regions 
during perihelion. 


3.2 Introduction 


The water production rate of comet 67P/Churyumov-Gerasimenko has been measured by 
a variety of instruments on board Rosetta. In June 2014, MIRO (Microwave Instrument 
for the Rosetta Orbiter; see Section 3.3 for more details) measured the water produc- 
tion rate at a heliocentric distance of 3.9 AU to be 107° molecules/second, increasing to 
4x10? molec/s by late August (Gulkis et al. 2015). Only Bockelee-Morvan et al. (2010) 
and de Val-Borro et al. (2014) have determined water outgassing rates from comets at 
larger heliocentric distances. Also in August 2014, Lee et al. (2015) demonstrated that 
the outgassing across the surface can vary by a factor of 30 from 0.1x10? molec/s/sr 
to 3.0x10? molec/s/sr. In the next month, MIRO estimated the production of the wa- 
ter isotopologue H5!6O as (4.9+2.5)x10°° molec/s on September 7 at 3.41 AU from the 
Sun (Biver et al. 2015). The ROSINA-COPS (Rosetta Orbiter Spectrometer for Ion and 
Neutral Analysis - COmetary Pressure Sensor) instrument found production rates in the 
range of 8.7x10? - 1.1x10?6 molec/s between August and November 2014, which then 
increased by a factor of 2 between November and January 2015 (Bieler et al. 2015). 
Fougere et al. (2016) present results from ROSINA-DFMS (Double Focusing Mass Spec- 
trometer), which show the outgassing increasing from «1076 molec/s at 3.5 AU (August 
2014) to >10?’ molec/s at 1.5 AU from the Sun (May 2015). A mean production rate 
of 8x10°° molec/s agrees with results from VIRTIS-H (Visible InfraRed Thermal Imag- 
ing Spectrometer) between November 2014 and January 2015 (Bockelée-Morvan et al. 
2015), and data in April 2015 from VIRTIS-M suggest that the rate increases to about 
10?’ molec/s (Migliorini et al. 2016). The ROSINA spectrometer is an in situ instrument 
unlike MIRO and VIRTIS, which are remote sounding instruments. 

Together, these data demonstrate the increasing activity of the comet towards peri- 
helion. In this work, we present the long-term evolution of the water production rate 
of comet 67P from MIRO observations as it approaches and recedes from the Sun, be- 
tween August 2014 and April 2016, with the aim to spatially resolve the water outgassing. 
Through this, we can analyse the behaviour of the comet with respect to heliocentric dis- 
tance and the change in activity across the surface. We use a lookup table method for the 
inversion of line areas into column densities and production rates and hence derive local 
effective Haser production rates of water (defined in Section 3.4.2) from the line area ra- 
tios of H5/5O and H5!5O. This approach is simple, producing results that are consistent 
with other measurements from Rosetta instruments along with ground-based observa- 
tions. Our methodology which enables large quantities of data to be inverted quickly and 
reliably, allowing the production rates from the 21-month MIRO dataset to be calculated 
in minutes and enabling us to study the outgassing evolution with good spatial resolution. 
With the spatially resolved data, we can study the processes driving the gas activity and 
investigate how the shape and illumination modulate this behaviour. The lookup table 
method is fully described in Section 3.4 along with a thermal sublimation model to help 
interpret the results, which are given in Section 3.5. The work is then summarised in 
Section 3.6. The next section gives a brief description of MIRO, the observations and the 
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shape model. 


3.3 MIRO observations 


The Microwave Instrument for the Rosetta Orbiter (MIRO) is a small millimetre/sub- 
millimetre spectrometer with a 30 cm offset parabolic reflector telescope and two hetero- 
dyne receivers. The millimetre receiver operates at a centre-band frequency of 188 GHz 
(1.6 mm wavelength) while the sub-millimetre receiver is tuned at 562 GHz (0.5 mm). 
Each frequency band contains a single broadband channel for the measurement of near 
surface temperatures. The continuum data have been previously analysed by Schloerb 
et al. (2015) and for polar night conditions by Choukroun et al. (2015). For our work, 
we are interested in the data from the Chirp Transform Spectrometer (CTS; Hartogh and 
Hartmann (1990)), which is connected to the sub-millimetre receiver. 

The CTS has 4096 channels with a spectral resolution of 44 kHz, corresponding to a 
resolution of v/Av = 10’ for the sub-mm channel. It can observe three isotopologues of 
water (H,!°O, H,!’0, and H,!8O), three methanol lines, carbon monoxide and ammonia at 
frequencies of 547 GHz to 580 GHz. A high spectral resolution is necessary for observing 
the narrow, Doppler broadened spectral lines of the low temperature, low pressure coma. 
MIRO collects spectral data every 30 seconds in a frequency switched mode. For the full 
technical details of the instrument, see Gulkis et al. (2007). We focus on the MIRO data 
collected from 3.62 AU in August 2014 to 1.24 AU at perihelion on August 13, 2015, and 
on the outbound leg to 2.02 AU in April 2016. In this time interval, the nucleus of the 
coma was fully resolved during the nadir observations. 

We use a digital shape model of the comet (SHAPS v1.2, Preusker et al. (2015)) 
with 200,000 triangular facets and the SPICE software provided by the Jet Propulsion 
Laboratory (JPL; Acton (1996)) to calculate the intersection geometry. With the shape 
model, the SPICE software and the beam pointing provided by the project, we can track 
the movement of the spacecraft and the comet to determine the facet-averaged position of 
the MIRO beam, solar illumination, solar zenith angle and viewing angle. 


3.4 Method 


We employed an efficient method using lookup tables to determine the water production 
rate of comet 67P. The basic idea of lookup table inversion is straight forward. We precal- 
culated a table of water production rate, Q, versus line area via a forward radiative transfer 
model (Section 3.4.2). Once the tables are calculated, the radiative transfer calculations 
are no longer required, and the measured line area ratios entered into a lookup table obtain 
the corresponding water production rate. In the application to the MIRO data, we used 
additional information to provide more unique mapping (using more than one table), such 
as the observed Doppler shift and continuum temperature. It is necessary to use more 
than one lookup table to test the sensitivity of the production rate measurements to each 
variable (Doppler shift and continuum temperature). We considered only nadir absorp- 
tion spectra in this work and treated the line area ratio, H590/H;!*O, as the fundamental 
variable. Below, the details of the methodology and measurement selection are provided. 
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Our approach is different from work by Biver et al. (2016), who have tracked the 
change in the global water production of 67P with heliocentric distance, as well as other 
molecules. In this work, we looked exclusively at nadir observations of the comet to 
derive local production rates and we neglected any geometrical effects. We then used the 
local production rates to track the evolution of the outgassing in each region as function 
of heliocentric distance. 


3.4.1 Water absorption lines 


In this analysis we focus only on the nadir absorption lines for the two water isotopes, 
H,'°O and H,'8O owing to their different opacities, and because they are routinely de- 
tected in the time period of interest (August 2014 to April 2016). Examples of these 
absorption lines can be seen in Gulkis et al. (2015), Biver et al. (2015), and Lee et al. 
(2015) (Figure 3), where the latter makes a distinction between the appearance of nadir, 
off-nadir, and limb viewing spectra. Processing the spectra involves several steps. First, 
we excluded spectra with low signal to noise (the line area must be greater than 2 K km/s 
for both isotopes), and then we excluded any emission lines, including limb and partial 
limb spectra. The 2 K km/s value for the noise limit was determined by calculating the 
average line area for continuum emission spectra without a line feature and with only 
random noise. The observations were averaged into 30-minute time intervals to reduce 
random noise; for example, in April 2016, the noise level decreased from 1 K km/s to 
0.14 K km/s as we increased the time interval for averaging the spectra from 30 seconds 
to 30 minutes. Subsequently, we smoothed the spectra with an 8 point box car filter. This 
smoothing has a minimal effect on the determination of the spectral line areas (<0.1%), 
but allows for a better determination of the Doppler shift by the automatic lookup method. 
The Doppler shift was simply calculated from the velocity shift of the minimum value of 
an absorption line. 

We found the line area by integrating in the region where the line forms, between -2 
and 1.1 km/s, which is wide enough to ensure that broadened lines are not cut off even for 
the wide H5!6O lines in the case of high gas production. The estimates for the line area 
error come from the random noise in the line wings (greater than 1.2 km/s) and calculated 
with the equation, 


ao = VN Av OT, (3.1) 


where N is the number of data points between -2 and 1.1 km/s, Av is the average veloc- 
ity interval between each point in the spectrum, and or is the standard deviation of the 
antenna temperature points in the wings (>1.2 km/s) 

In Figure 3.1, the time evolution of the nadir line areas for the two isotopes are pre- 
sented. For both isotopologues, the activity increases towards perihelion owing to the 
increase in outgassing, and decreases after perihelion, but with a delay (discussed later). 
The line area uncertainties due to random noise are small, typically 0.2 K km/s for both 
isotopologues with measurements in the range of 13 - 385 K km/s for H,!6O and from 2 
-135 K km/s for H,!?O (see Figure 3.1). 

The line area ratio is a sensitive tracer of opacity and column density in the coma 
(Figure 3.1). Providing that the optical depth is less than unity, the line area is linearly 
related to the column density. However, as the column density increases further, line area 


50 


3.4 Method 





400 























j: 


T na 
$ 


| IR TEN IE ‘ 
Oe UH Maa La dr | 
5 


9 Sep 14 Nov 14 Jan 15 Mar 1 May 15 Jul 15 Sep 15 Nov 15 Jan 16 Mar 16 May 16 











Figure 3.1: Time evolution of the line area of H5!5O (top panel), H5!*O (middle panel), 
and the H5'^0/H5?O line area ratio (bottom panel) from August 2014 to April 2016. 
The thick dashed line represents perihelion, which occurred on August 13, 2015, and the 
thin dashed lines represent when 67P was 1 AU from closest approach, pre- and post- 
perihelion. 


growth loses this linearity and the dependence becomes Vin N, where N is the number 
density (Goldsmith and Langer 1999). The line is said to be saturated, as the amplitude 
does not change with increasing column density. The H,!°O and H,!8O lines observed by 
MIRO follow this pattern. Because these two lines have different opacities, the H,!°O line 
is always saturated with an opacity >10 (Gulkis et al. 2015) whereas the H5!*O line starts 
out with a small amplitude and is optically thin (opacity «0.1) but then grows with the 
column density. By combining them in a ratio we produce a more favourable line growth 
curve with greater sensitivity as we approach the saturated part of the growth curve. 

In Figure 3.1, the ratio indicates that in August 2014 the H,!°O is optically thick as the 
highest values are around 60-70, much lower than the isotopic ratio of 498.7 (terrestrial 
value). The ratios further decrease towards the perihelion as even the H2!*O line becomes 
optically thick, and the ratio reaches values as low as 2.5. About 90 days after perihelion 
the ratio grows again significantly as the water production goes down. 


3.4.2 Generating the lookup tables 


We performed the lookup table calculations with a non-LTE forward model. For radia- 
tive transfer in comets, two different computational methods are typically used: a locally 
approximated escape probability approach (Bockelée-Morvan 1987, Litvak and Kuiper 
1982) and a direct ray transfer Monte Carlo approach (Hogerheijde and van der Tak 2000). 
Both methods have been found to produce similar results numerically (Zakharov et al. 
2007) and in comparison to observations (Hartogh et al. 2010b, de Val-Borro et al. 2010). 
Building a lookup table involves hundreds of calculations, which with the Monte Carlo 
approach, would take weeks to compute. We instead used the escape probaility approach 
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following Lee et al. (2015) (see references therein for detailed description and accuracy), 
which required substantially less computational time. 

The atmospheric inputs are defined radially over 500 km from the nucleus surface 
with grid points (integer index 1) calculated at 


R4C! PRENT (3.2) 

2 
where R, is the radius of the comet, taken as 2 km, multiplied by a scaling constant, 
C=1.02, to define 278 logarithmic grid points over 500 km. The temperature profile fol- 
lows an adiabatic expansion. The inputs are the surface temperature parameter, T,, and the 
terminal temperature, which is estimated from the line cores of the H,!°O lines observed 
by MIRO and fixed at 30 K. The velocity profile is modelled by a hyperbolic tangent func- 
tion and increases from zero to the terminal expansion velocity, v,,, , as described in Lee 
et al. (2015). The radial molecular number density profile is dictated by the Haser model 
(Haser 1957) which depends on the production rate, Q, expansion velocity profile and 8 
parameter, which is also taken from Lee et al. (2015). The 6 parameter is the photodis- 
sociation rate with a value of 1.4x10~> s~! at a distance of 1 AU, but it is not critical for 
deriving Q. In addition, we assume the H,!°O/H,'?O ratio and the ortho-para ratio to be 
500 and 3, respectively, in our forward calculations, as in previous works. In total, three 
atmospheric parameters are varied to generate two lookup tables for the water production 
rate: the surface temperature, T, (K); the terminal expansion velocity, v.., (km/s); and Q 
(molec/s). 

In the first table, the production rate varies from 10? to 10?? molecules/s over seven 
different velocities (0.4 - 1 km/s) and produces the curves given in Figure 3.2. To avoid 
interpolation within the tables, a ninth order polynomial is fitted through each curve with 
an accuracy of >0.02%, which then provides an analytic expression for each function. 
The errors on the ratios are treated as upper and lower bounds which are also calculated 
with the same polynomial equation to give the errors on the production rates. Only the 
production rate and velocity parameters are varied; the surface temperature parameter, T,, 
is kept at 200 K in the temperature profile. 

In order to correct the measurements for temperature sensitivity, we created a second 
lookup table. It is made in a similar way but now T; is allowed to vary from 50 - 300 K, 
while the velocity parameter is held at a constant value of 0.7 km/s. Variations in T, would 
change the line contrast. The production rate ranges from 10% to 10? molecules/s. Here, 
for each value of the production rate, the temperature is increased and the ratio between 
the two absorption lines is calculated again. Figure 3.3 shows these functions with line 
area ratio versus continuum temperature. Each contour line represents a different produc- 
tion rate with the lowest rate (1074 molecules/s) at the top of the figure and the highest 
rate (107? molecules/s) at the bottom. We took the measured MIRO data, overplotted 
them onto these curves, and performed a linear interpolation between the lines to find the 
water production rate. This was also performed for the line area ratio errors as before to 
produce upper and lower bounds on the production rate estimates. For each production 
rate, as the continuum temperature increases, the ratio between the two absorption lines 
also increases. 

The weaker H5!5O line is formed mainly near the nucleus. The temperature decreases 
from the surface, so there is a smaller contrast between the H5'?O line core and the con- 
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Figure 3.2: Log-log plot of the H50/H5!*O absorption line area ratio as a function of 
production rate for 7 different terminal expansion velocities (0.4 - 1 km/s). For each 
velocity, as the production rate is increased, the ratio between the two absorption lines 
decreases from about 500 to approximately 1. The measured ratio and expansion velocity 
from the MIRO dataset can be used to derive estimates for the production rate from these 
curves. 


tinuum temperature. On the other hand, the opaque H5!9O line forms further from the 
nucleus, at altitudes where the temperature profile is approaching the terminal tempera- 
ture (estimated to be 30 K, as discussed earlier in this section), and therefore has a larger 
contrast. Increasing the continuum temperature therefore increases the ratio of line areas 
between the strong and weak lines. 

A second lookup table is needed to test the sensitivity of the production rates to tem- 
perature, which was assumed to be constant in the first lookup table. The effect of this 
second table is a small correction which typically increases the production rates derived 
from the first lookup table by about 20-2596. The relative uncertainty on the production 
rates with these lookup tables are in the range of 30-5096. This uncertainty grows though 
with increasing Q when Q>10°. This particular method is not suitable for deriving water 
production rates where Q»10? as the line area versus production rates loses sensitivity 
above this threshold because both the H,!°O and H5!*O lines become saturated. As illus- 
trated in Figures 3.2 and 3.3, the close spacing of the contour lines limits the accuracy on 
determining Q from the lookup tables, but as we want to examine the relative behaviour 
of the production rate from the comet, this method still provides fast and reliable infor- 
mation. 

The Haser model gives the global production rate for a spherical comet but our cal- 
culated results are local production rates determined for the viewing position of the beam 
on the comets surface. Each measurement is then a local effective Haser production rate, 
meaning that the localised beam viewing spot appears to be outgassing as if it was a spher- 
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Figure 3.3: Log-linear plot of the H590/H5!*O absorption line area ratio as a function of 
continuum temperature for 101 different production rates (10?^ - 10? molecules/s from 
top to bottom). For each production rate, as the continuum temperature increases, the 
ratio between the 2 absorption lines also increases. The measured ratio and continuum 
temperature from the MIRO dataset (grey points) can be used to derive estimates for the 
production rate from these curves. 


ical comet with a single global production rate. This local effective Haser production rate 
is proportional to the derived MIRO column density. 


3.4.3 Thermal sublimation model 


For later interpretation of results, we now describe a basic thermal model for predicting 
the water outgassing rate, which relies on the comet A model in Keller et al. (2015). We 
used the SHAP2 shape model (Sierks et al. 2015), which has fewer facets than SHAPS5. 
We performed the described calculations for each facet of the digital shape model, except 
that we neglected self-heating mechanisms and shadowing. Keller et al. (2015) com- 
pare their A, B, and C models with and without self-illumination and show that self- 
illumination mainly plays a role when the comet is far away from the Sun (23.5 AU). 
In the time period examined here, the difference in water production rate is negligible if 
self-illumination is ignored. The sublimation rate can be expressed as 
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(1-A,)I = eoT* + Z(T)Lice, (3.3) 


where A, is the bolometric Bond albedo, taken as 0.01, I is the solar flux received by a 
facet, € = 0.9 is the emissivity, o is the Stefan-Boltzmann constant, T is the facet tem- 
perature, Z is the facet sublimation rate, and L;.. = 2.6x10° J/kg is the latent heat of 
sublimation for water ice, which we take to be constant. The sublimation rate, water 
vapour pressure, and thermal velocity are given by (Equation 5 from Panale and Salvail 
(1984)) 


Z(T) 2 2P(T)/(va) (3.4) 
P(T) 2 3.56 x 10? exp(-6141.667/T) [kgm ' s ?] (3.5) 


va (T) = NSRT nu (3.6) 


The surface temperature in this model plateaus at around 200 K, which is the same 
result found for Model A in Figure 4 of Keller et al. (2015). It is worth noting the non- 
linear temperature sensitivity of the water vapour pressure and hence the sublimation 
rate. Increasing the temperature from 140 K to 200 K also increases the vapour pressure 
from 10°’Pa to 10°!Pa so that even small variations within this range of temperatures has 
profound effects on the sublimation rate, 

The production rate is proportional to the sublimation rate: 


Q - AZ (3.7) 


where A is the area of the active surface; we assume all facets to be active. By investigat- 
ing how and why the sublimation rate changes in this thermal model, we can infer how 
the production rate changes as well. 

The shape model gives the average daily illumination of every facet and Equation 3.3 
allows for the calculation of the temperature and sublimation rate for each facet. With 
this information, the average temperature and sublimation rate of the facets across the 
surface can hence be determined for a given day with these averages weighted by the area 
of each facet. Using the same latitude bins described in Section 3.5.2, we also calculate 
the average sublimation rate for six subset regions. It is possible to go one step further, 
and divide the grid into 20 longitude bins and hence calculate the sublimation rate of 120 
sections across the surface. The model is used from 250 days before perihelion to 250 
days afterwards to track the behaviour of the sublimation rate. The sublimation rate at 
each point in time is plotted to determine the r, exponents before and after perihelion. 
The dataset from the thermal model is trimmed so that we only calculate the exponents 
over the period of time for where there are observations. 


3.5 Results 


3.5.1 Water production rate 


Figure 3.4 shows the derived local effective Haser production rates of water from 250 days 
before and after perihelion. For comparison, we overplotted other measurements made by 
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instruments on board Rosetta (red) and ground-based observations (blue). In addition, 
modelled production rates from Keller et al. (2015) are shown for completeness. There 
is a general qualitative agreement between all the data included; the MIRO estimates 
derived here from the lookup tables are consistent with early results (-400 to -100 days 
from perihelion) from VIRTIS-H (Bockelée-Morvan et al. 2015), VIRTIS-M (Migliorini 
et al. 2016, Fink et al. 2016), ROSINA-DFMS (Fougere et al. 2016), and ROSINA-COPS 
(Bieler et al. 2015) along with previous MIRO estimates (Gulkis et al. 2015). Around 
120 days before and after perihelion, the derived MIRO estimates and ground-based ob- 
servations (Bertaux et al. 2014, Crovisier et al. 2002, Feldman et al. 2004, Hanner et al. 
1985, Ootsubo et al. 2012, Schleicher 2006) are in rather good agreement, although these 
were made for previous perihelion passages of comet 67P. In addition, a power law line 
shown in Fougere et al. (2016) and a model result for a comet with a perpendicular spin 
axis and an active area of 296 from Cowan and A'Hearn (1979) cross through the derived 
MIRO production rates. Hansen et al. (2016) (not shown in Figure 3.4) have derived the 
long-term variation in the global water production rate of comet 67P from ROSINA data, 
but our local production rates are generally lower than their findings, as are the VIRTIS 
measurements (Bockelée-Morvan et al. 2015). We find that the largest local production 
rate for a single observation is (1.42+0.51)x1078 molecules/s on August 29, 2015, and 
this is a factor of 3 lower than the maximum outgassing given by Hansen et al. (2016). 
The difference may be because we derive local Haser production rates whereas their work 
calculates global production rates. Additionally, the two instruments use different acqui- 
sition techniques: ROSINA is an in situ instrument whereas MIRO measures line-of-sight 
spectra. Hansen et al. (2016) discuss how their result is a factor of two larger than results 
from VIRTIS, which is also a line-of-sight instrument. 

The scatter seen in the local water production rates are indicative of the variation in 
latitude from the subsolar point of each observation. In Figure 3.4, the peaks in the local 
effective Haser production rates can be seen to occur when MIRO views latitudes close 
to the subsolar latitude and conversely, the dips correspond to times when latitudes away 
from the subsolar latitude are observed. 

A 14-day moving average smoothing is applied to the MIRO data to give a qualitative 
description of the mean behaviour over the observational period, shown in orange and 
in Table 3.1. The moving average assumes that day- and night side observations are 
equally weighted. The water outgassing increases from (3.16+1.00)x10*° molecules/s 
in the first averaging bin (i.e. 231 days before perihelion, ~2.4 AU from the Sun) to 
(5.05+2.46)x 10?’ molecules/s a week before perihelion. Then 245 days after perihelion 
(the last averaging bin), the water production rate drops to (1.93+0.46)x 107° molecules/s. 
The bias in the latitudinal sampling means that the moving average is close to, but does 
not represent, the global production rate: the average may over- or underestimate the true 
global value depending upon the distribution of the observed latitudes. Nevertheless, the 
qualitative agreement between the moving average and global production rates from other 
instruments and ground-based observations is good, and these global measurements also 
fall within the point cloud of the local effective Haser production rates derived in this 
work. 

The black lines for models A and D represent a comet with an even, homogeneous 
covering of ice across the surface, whereas models E and F represent a comet with a 
spotted, heterogeneous distribution of ice (Keller et al. 2015). Model A reproduces the 
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Figure 3.4: Top panel: Calculated local effective Haser production rates of water from 
250 days before and after perihelion from the MIRO dataset (grey). Results from other 
Rosetta instruments of global production rates (red) and ground-based observations from 
previous perihelion passages (blue) are shown for comparison as well as results from 
modelling (black lines). A two-week moving average of the MIRO data is also included 
(orange line). Bottom panel: The latitude for each MIRO observation (grey points) and 
the position of the subsolar latitude over time (black line). 


production rates reasonably well at perihelion but the model values are higher in compari- 
son to other models during the approach to perihelion as also noted in Keller et al. (2015), 
who found that results were enhanced by a factor of 16. The case is the same for Model D, 
which additionally considers thermal inertia. In these two models, all facets are deemed 
to be active. For models E and F, which have different values of thermal inertia, some 
facets are switched off and the abundance of ice is increased elsewhere. The production 
rate around perihelion for these models does not change in comparison to models A and 
D, but the outgassing is reduced pre and post-perihelion. This is more consistent with 
the derived MIRO results. Hence, qualitatively, the average two week production rate 
(orange line) appears to be more consistent with the heterogenous surface outgassing as 
described in models E and F in Keller et al. (2015). The authors of this paper also found 
that these two models were more consistent with pre- and post-perihelion ground-based 
observations (Schleicher 2006, Hanner et al. 1985, Ootsubo et al. 2012) than the other 
models. 


The data also clearly demonstrate that the maximum production is shifted away from 
perihelion. To investigate this, linear fits are made through the pre- and post-perihelion 
data that converge on a solution where the two lines intercept. In order to do this, it is 
necessary to restrict the dataset, as the H2!*O line is small at distances greater than about 
2.5 AU (before February 2015). As a consequence, the H,'?O line area can be below 
the desired signal-to-noise level and without the lowest line areas, we sometimes miss 
the highest ratios and therefore the lowest water production rate measurements. This can 
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be seen in Figure 3.4 in the apparent plateauing of the MIRO results 180 days before 
perihelion (orange line). For the offset calculation, we hence restrict the dataset to look 
at results within 1 AU of perihelion (March 2015 - February 2016), where this effect is 
not prevalent. The intercept from linear fitting in log Q-log r; space gives the solution as 
34+10 days after perihelion. Previous work by Bondarenko and Medvedev (2011) and 
Ferrin (2010) found from cometary light curves that the offset occurred approximately 
33 days after perihelion. Schleicher (2006) also found an offset in the peak outgassing 
from perihelion by 1 month using OH gas production rates. Furthermore, observations 
of 67P around its three previous perihelion visits (1996, 2002 and 2009) suggest that 
the activity peaks 16 + 5 days after perihelion and then again 40 days after perihelion 
(Bertaux et al. 2014). The offset found here is consistent with Bondarenko and Medvedev 
(2011), Ferrin (2010) and the second peak from Bertaux et al. (2014), but we are unable to 
resolve the first peak in this latter work. There is a data gap in mid-September 2015 when 
the spacecraft goes above 500 km, which limits the determination of the offset. Above 
500 km, the MIRO beam is larger than the comet and includes contributions from the 
background and the nucleus, so MIRO does not collect pure nadir absorption spectra in 
this period. 


3.5.2. Behaviour with heliocentric distance 


Linear fitting in log Q-log r, space around perihelion of the local effective Haser produc- 
tion rates, allows for the estimation of the change in water outgassing with heliocentric 
distance by evaluating the r, exponents, which we define as a. We find that these have 
values of -3.8 + 0.2 and -4.3 + 0.2, within 1 AU pre- and post-perihelion, respectively, 
taking into account the offset described in the last section. This means that the local 
effective Haser production rate, Q, across the entire surface is roughly symmetrical and 
proportional to the heliocentric distance, r;,, as approximately Q « r,*. The thermal mod- 
elling described in Section 3.4.3 calculates a for pre- and post-perihelion as -3.5 + 0.01 
and -3.9 + 0.03 for the entire surface. 

We further investigated the behaviour by calculating a at different latitudes. We split 
the dataset into six latitude bins: 90° to 25°; 25° to 0°; 0° to -20°; -20° to -30°; -30° to 
-45°; and -45° to -90°. These bins were chosen so that each bin has the same number of 
data points (200). As Figure 3.5 shows, a has a latitudinal dependence: the slopes steepen 
from north to south for pre- and post- perihelion, from r, exponents of about -2 to -6 (right 
panel). The error bars reflect the standard error of the linear fitting. In addition, the data 
indicate that there is a line intersection where the post-perihelion slopes are steeper than 
the pre-perihelion slopes for latitudes greater than 0°, but at more southern latitudes, the 
slopes are steeper during pre-perihelion than post-perihelion. 

To explain this behaviour, we turn to the results of the thermal sublimation model 
(Keller et al. 2015) which are shown in the left panel of Figure 3.5. This panel shows 
the @ values for the sublimation rate before and after perihelion in the six latitude bins. 
There is also a steepening of the slopes from north to south similar to the trend seen in 
calculated production rates. The standard deviation reflects the variation in longitude for 
each latitude bin. 

The steepening of the slopes from north to south in the MIRO data and the model can 
be explained by the rotation axis of the comet which has been found to be tilted by about 
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Date Days from | Average local production 
perihelion rate (x107’molec/s) 
25/12/2014 -231 0.32+0.18 
08/01/2015 -217 0.37+0.15 
05/02/2015 -189 0.35+0.13 
19/02/2015 -175 0.37+0.11 
05/03/2015 -161 0.68+0.43 
19/03/2015 -147 0.86+0.56 
02/04/2015 -133 0.57+0.32 
16/04/2015 -119 1.03+0.53 
30/04/2015 -105 1.02+0.51 
14/05/2015 -91 1.55+0.74 
28/05/2015 -77 1.97+0.71 
11/06/2015 -63 2.29+1.00 
25/06/2015 -49 2.42+1.34 
09/07/2015 -35 3.16+1.60 
23/07/2015 -21 4.07+2.10 
06/08/2015 -7 5.05+2.46 
20/08/2015 7 3.68+2.14 
03/09/2015 21 6.11+2.66 
17/09/2015 35 3.81+2.06 
01/10/2015 49 6.16+1.91 
15/10/2015 63 3.98+1.99 
29/10/2015 77 4.824152 
12/11/2015 91 2.10+1.12 
26/11/2015 105 3.03+1.37 
10/12/2015 119 1.98+0.96 
24/12/2015 133 0.93+0.75 
07/01/2016 147 1.25+0.55 
21/01/2016 161 0.99+0.51 
04/02/2016 175 0.77+0.27 
18/02/2016 189 0.57+0.25 
03/03/2016 203 0.22+0.06 
17/03/2016 217 0.21+0.09 
14/04/2016 245 0.17+0.05 

















Table 3.1: Fourteen-day moving average applied to the local effective Haser production 
rates of water calculated from the MIRO data. 
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50° to the orbital axis and so the subsolar latitude moves from approximately 40°N in 
August 2014 to -50°S in September 2015, just after perihelion (Keller et al. 2015). The 
northern latitudes are well illuminated when the comet is further away from the Sun (3.5 
AU) but not illuminated during perihelion when the received solar flux increases. As a 
result, the surface temperature of the northern hemisphere does not change that much in 
the period for which there are observations and, therefore, neither does the sublimation 
rate, which is very sensitive to changes in temperature, as already discussed. Conversely, 
the southern latitudes are not illuminated when the comet is more than 1.8 AU from the 
Sun, but very well illuminated at perihelion. The change in the surface temperature of 
the southern hemisphere is quick and as a result, so is the sublimation rate. The resulting 
exponents for the southern hemisphere are therefore less than in the northern hemisphere. 

The line intersection described above is also present in the model, as the post-perihelion 
slopes are steeper than the pre-perihelion slopes at northern latitudes but the situation 
is reversed at southern latitudes. In the southern hemisphere, during the pre-perihelion 
phase, the surface temperature increases very suddenly before perihelion but the decrease 
in temperature afterwards is slower, so the slope for the sublimation rate is steeper before 
perihelion than afterwards. For the northern hemisphere, during the approach to peri- 
helion, the sublimation rate changes very little and only increases slightly, whereas the 
decrease is more pronounced after the approach to perihelion. 

The data and model results show that o varies with latitude and that a single a value 
representing the whole comet surface may not be of fundamental significance. The pro- 
duction rate is a strong non-linear function of temperature, which in turn, depends upon 
the shape, morphology, and obliquity of the comet, along with the rate of change of illu- 
mination; therefore a single a value on its own may not reveal much about comet activity. 
It appears that the same distribution of ices over a sphere would produce completely dif- 
ferent results for a. The latitudinal dependence is more important and it can help to reveal 
areas of the comet that are particularly active or inactive. This can only be investigated 
by in situ measurements as collected by MIRO. 

There are some differences between the a values from the measured production rates 
and the calculated sublimation rates. However, all of the measurements from the MIRO 
data are within 1 standard deviation from the model results except for the post-perihelion 
points north of -20 latitude and the pre-perihelion point north of 25°, which are within 
2 standard deviations and have overlapping error bars. For most of points, the standard 
deviation is between 0.7 and 1.5 but the more southern latitudes have a much larger vari- 
ation in longitude pre-perihelion. We have already stated that the more southern latitudes 
experience a sudden increase in temperature but owing to the local morphology of the 
surface it can be even more extreme at some longitudes than at others. There is then, a 
broad range of possible œ values across the longitudinal plane. Despite the simplicity of 
the thermal sublimation model that we used here, it captures and reproduces much of the 
behaviour seen in the MIRO data to a good degree of accuracy. One of our modelling 
assumptions is that the active area is constant across the comet, which is probably not 
true (Keller et al. 2015). In addition, the highest (770^) and lowest latitudes (<-80°) are 
not really sampled by MIRO and could be excluded from the model; also there may be a 
longitudinal sampling bias in the data, although this is accounted for in the error bars of 
the model. 

It should be noted that there is a slight difference between how the sublimation rate 
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Figure 3.5: (Left) Modelled pre- and post-perihelion slopes (o) of water sublimation rate 
as a function of latitude averaged over several bins (see text for detailed description). 
(Right) Pre- and post-perihelion slopes of water production rate from the data derived 
from MIRO measurements. The error bars on the modelled data reflect longitudinal vari- 
ability, which is very large in southern hemisphere owing to rapid temperature variations. 
The MIRO data error bars also reflect the data variability at 1o level. 


and the MIRO production rates are calculated. The sublimation rate is calculated from 
the thermal model, which takes into account a whole latitude band, spanning the day- 
and night sides, and all longitudes. On the other hand, the production rates only come 
from a spot on the nucleus where the MIRO beam is observing, which could be day or 
night, and only at a certain longitude. However, with an appropriate level of sampling 
(approximately 200 MIRO measurements per latitude bin) enough of the surface is ob- 
served to average out the effects of longitudinal variance. A more sophisticated approach 
where only the facets of the shape model observed by the beam are used to calculate the 
sublimation rate would be more representative of the actual measurements, but the sim- 
ple thermal model has captured the essence of the results and can explain some of the 
features. 


3.5.3 Regional variations 


The large dataset obtained through the fast inversion can be also spatially partitioned to 
study regional activity variations using the latitude and longitude co-ordinates of the beam 
position. There are 26 identified regions on 67P defined by El-Maarry et al. (2016) and 
we can identify the regional source of the produced water from the shape model using the 
average beam position for each measurement. 

First, however, if we assume that molecules move radially away from the comet, then 
we must restrict our dataset to measurements with a low viewing angle (i.e. the angle 
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between surface normal and beam direction) to ensure that we are only observing the pro- 
duction from the region that we are interested in. Using a large viewing angle means that 
the MIRO beam cuts through layers of the coma where the water production originates 
from other areas. We hence choose measurements with a viewing angle of less than 30° 
to reduce the uncertainty in the regional location of the MIRO beam. This is calculated 
from the SHAPS5 shape model and it is the average viewing angle of all the facets within 
the beam. We assume that all the contributions to the local effective Haser production rate 
in the line of sight originate from the MIRO footprint beam. This assumption is only an 
approximation as it is possible that a particular line-of-sight column density is influenced 
by other regions owing to the complex nucleus shape and larger phase angles. These 
potential contributions should be carefully quantified, but it is out of the scope of this 
work. 

The results are given in Figure 3.6, showing regional variations in the water outgassing 
across the surface of the comet. The average daily solar flux for each region is given in the 
lower panel for comparison. It is calculated from the solar flux constant and the average 
cosine solar zenith angle for each facet within the region. 

For all regions, the water production rate increases as 67P approaches perihelion and 
decreases afterwards. The change in production rate is not uniform though. The southern- 
most regions (Anhur, Geb, Bes, Sobek, Neith and Wosret) experience a dramatic change 
in outgassing and are the origins of the largest measured production rates around peri- 
helion. As already mentioned, the subsolar latitude moves south and goes below -50° 
latitude just after perihelion. These southernmost regions are therefore the most illumi- 
nated and receive the largest amount of flux at perihelion (between 400 and 600 W/m?). 
This may result in a stronger depletion of water and potentially other volatiles in these 
regions. 

There are four regions (Atum, Anubis, Khonsu, and Imhotep) at about -30° latitude 
that are also reasonably well illuminated around perihelion (with fluxes between 200 
and 400 W/m?) and have reasonable local effective Haser production rates (as high as 
~ 9x10?’ molec/s), although these are slightly lower than most southern regions. As we 
look further north, to the equator, the water production rates from Ap, Anuket, Maftet, 
Bastet, Khepry, and Aker are considerably smaller. Due to the obliquity, the northern re- 
gions spend most of their time facing away from the Sun at perihelion (whilst the southern 
hemisphere is well illuminated) so the received flux is lower. This explains why the HO 
outgassing is lower as we look towards northern latitudes. Indeed, for the regions around 
the north pole (Seth, Ash, Aten, Babi, and Hapi) as well as on the head lobe (Hathor, 
Hatmehit, Ma’at, Nut, and Serget), the production rates are found to be about 50% lower 
than in the south at perihelion. 

Generally, the regions with high received solar fluxes (2200 W/m?) have consistently 
high water production rates (see two leftmost panels in Figure 3.6) and that insolation 
explains most of the activity. In the other regions, the flux is lower («200 W/m?) and the 
production rates are suppressed. It is notable that, for all regions, the largest production 
rates are shifted towards perihelion even though in some areas, the received flux is low: for 
example, in Seth and Ash, the flux drops below 100 W/m? but several of the production 
rate measurements are as high as those seen in the well-illuminated southern regions. 
Possibly in regions like these, solar insolation is not the only driver of activity but it is 
hard to draw too many concrete conclusions. It is also possible that these large production 
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Figure 3.6: Change in water production rate 250 days before and after perihelion sorted 
by regions, as defined by El-Maarry et al. (2016). The first panel shows the results from 
the six southernmost regions, Anhur, Geb, Bes, Sobek, Neith, and Wosret. The largest 
production rates are found here, specifically from the Bes, Wosret, and Neith regions. The 
second panel contains the results from four regions located around -30° latitude: Atum, 
Anubis, Khonsu, and Imhotep. These show reasonable outgassing rates (as high as ~ 
9x10? molec/s) and Imhotep in particular is well sampled in this time frame. Six regions 
that straddle the equator are in the third panel: Ap, Anuket, Maftet, Bastet, Khepry, and 
Aker. These is an apparent drop here in the production rate around perihelion as we start 
to look at more northern regions. Seth, Ash, Aten, Babi, and Hapi - the five northernmost 
regions - are in the fourth panel. Even though these regions do not receive much solar flux 
due to the rotation axis of the comet, Seth and Ash are still quite active. The remaining 
five regions, Hathor, Hatmehit, Ma’ at, Nut, and Serqet are all located on the head lobe and 
shown in the fifth panel. The average daily solar flux for each region is shown underneath 


for comparison. 


rates seen in Seth and Ash include contributions in the line of sight of the MIRO beam 
from other outgassing areas. A more in-depth study of these and other measurements 
is needed to properly identify any particularly active or inactive areas from the MIRO 
dataset. 


3.6 Conclusions 


The line area lookup tables method is a fast and computationally inexpensive way to cal- 
culate the production rate of water from comet 67P/Churyumov-Gerasimenko. Our results 
are consistent with measurements from other Rosetta instruments as well as ground-based 
observations and produces the time dependence accurately. The uncertainty on each wa- 
ter production rate measurement is between 30% - 50%. Systematic errors dominate the 
uncertainty, as can be seen in Figure 3.2, for example, where the ratio curves converge to 
unity and are flat for large water production rates, so errors in the ratio can give large error 
bounds on these estimates. In fact, when Q > 107°, this method does not provide accurate 
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results as the uncertainty can be a factor of 2-3. The errors in the line area measurements 
from random noise are small, as detailed in Section 3.4.1. The retrieved data variabil- 
ity in a given time period is of the order of o = 2, which reflects true data variability 
due to illumination and nucleus inhomogeneities, but also includes modelling errors from 
simplifications. These however, do not contribute more than 50%. The rapid inversion 
approach enables us to obtain nearly complete data coverage on the nucleus over time, 
and study in detail the spatial distribution of outgassing activity and the evolution with 
heliocentric distance. In combination with the shape model, the dataset can be partitioned 
by location on the nucleus surface with latitude and longitude co-ordinates. With this, we 
can determine active regions of the comet and investigate why the activity changes. 
The following conclusions can be drawn from the results: 


e During the approach to perihelion, the derived results for the local effective Haser 
production rates of water are consistent with previous global production rates mea- 
sured by MIRO, ROSINA and VIRTIS. 


e The MIRO 14-day mean average changes by an order of magnitude during the 
500-day window around perihelion: 231 days beforehand, the outgassing rate is 
(3.2441.79)x 10° molec/s, a week before perihelion it is (5.05+2.46)x 107’ molec/s, 
and then 245 days after perihelion it decreases to (1.72+0.54)x10”° molec/s. 


e By integrating this average outgassing curve, the comet is found to have lost (2.4 + 
1.1)x10? kg of water in this time period, which for a dust-to-gas ratio y =4 (Rotundi 
et al. 2015) is (1.2 + 0.6)x10"° kg of material or 0.12 + 0.06 % of the comets mass, 
assuming a total mass of 1.0x 10P kg (Sierks et al. 2015). For comparison, Hansen 
et al. (2016) give the mass loss of water as 6.4x10?kg from 3.6 AU inbound to 3.0 
AU outbound and a mass loss of 4 - 8 x 10'° kg in this period. Both of these results 
are a little higher than what we report here. Additionally, Bertaux (2015) derived a 
water mass loss of (2.7+40.4)x10° kg per orbit from SWAN/SOHO measurements 
and Keller et al. (2015) calculated the comet would lose 6.5x10!? kg of water during 
an orbit within their Model A. This latter result is significantly higher than what is 
derived here for the MIRO results but this is to be expected: Model A in their work 
assumes a homogeneously active surface, which leads to production rates that are 
as much as 16 times higher than observations during the early perihelion passage 
and therefore leads to a higher estimate for the water production over one orbit. 


e The maximum local effective Haser production rate of a single observation occurs 
on August 29, 2015 with a value of (1.42+0.51)x1078 molec/s or (426 + 153) kg/s. 
The outgassing appears to peak 34+10 days after perihelion. These results are also 
consistent with ground-based observations of 67P from previous apparitions. 


e Linear fitting in log Q-log r; space finds that the change in the water production rate 
across the surface of the comet follows power laws of Q c r;?5*?? pre-perihelion 
and Q « r,**°? post-perihelion. A thermal sublimation model derives similar 
values for the exponents with values of -3.5 + 0.01 and -3.9 + 0.03, respectively. 
The MIRO results also indicate that there is a latitudinal dependence on the produc- 
tion rate for which northern and southern regions have r, exponents (œ), as seen in 


Figure 3.5. This can be explained by the orientation and orbit of the comet. The 
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sublimation rate model driven by the solar illumination indicates that the northern 
hemisphere does not experience dramatic illumination changes during perihelion 
approach and it is in shadow 20 days before perihelion. On the other hand, the 
southern hemisphere suddenly becomes illuminated at perihelion, causing a dra- 
matic increase in the sublimation rate. The subsolar latitude moves south quickly 
in the first half of 2015 but the return north is slower after perihelion. 


e The retrieved local effective Haser production rates can also be used to identify the 
most active regions on the comet, particularly those that experience an increase in 
illumination around perihelion. The southernmost regions, including Wosret, Neith, 
and Bes, are the most active regions in this period. Despite the fact that the received 
solar flux in the northernmost regions and on the head lobe drops significantly (see 
rightmost panels in Figure 3.6), the local effective Haser water production rates 
can still be relatively high; for example, in the Seth region just after perihelion, the 
production rate was measured to be ~ 9x 10?’ molec/s even though the average daily 
flux was ~ 50 W/m?. Further work is needed to quantify the effects of geometry 
and phase angle on the water column determination from these particular regions. 


As noted by Maquet (2015) and Guzzo and Lega (2015), comet 67P has had a chaotic 
orbital history. It was only brought onto an orbit with a perihelion distance of 1.3 AU after 
a close encounter with Jupiter in 1959. The comet has now experienced nine perihelion 
passages at this distance. Before the 1959 encounter with Jupiter, the estimated closest 
approach to the Sun was 2.7 AU. Even though this comet might have spent most of its 
time relatively far from the Sun, it is likely that the major observed geomorphological 
forms, the apparent layering, and the variations in lateral volatile ice are all a result of 
evolutionary processing (Guilbert-Lepoutre et al. 2016). In addition, the loss of volatiles 
is currently modulated by the actual nucleus shape and spin axis orientation (Keller et al. 
2015). Modern models of the long-term evolution of cometary bodies demonstrate that 
even a non-uniform surface albedo distribution produces varying thermal gradients, which 
leads to inhomogeneous ice distribution and composition (Guilbert-Lepoutre and Jewitt 
2011). Rosenberg and Prialnik (2010) also show strong evidence that an initial heteroge- 
neous distribution of internal porous ice patches strongly affects the production rates and 
mass loss. In this connection, it is conceivable that the chaotic path of 67P may have in- 
duced the radial and lateral distribution of the cometary refractory and volatile content in 
the comet, which could explain the observed activity pattern of different volatiles. Hueb- 
ner et al. (2006) provide a comprehensive overview of the effects of multi-stage injection 
into the inner solar system on the internal structure of the nucleus. Hence, establishing 
a possible relationship between the present state and the original state of the comet is an 
important but difficult and unresolved task, which includes the estimates of the fraction 
of volatiles that the comet 67P has lost so far. Appropriate modelling efforts are needed 
to yield insights into these questions. 
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heliocentric water production rates 
of comets 


This chapter has been accepted for publication by A&A as the journal article: Marshall, 
D., et al. "On the interpretation of heliocentric water production rates of comets." 


4.1 Summary 


We investigate the influence of three basic factors on the water production rate as a func- 
tion of heliocentric distance: nucleus shape, the spin axis orientation, and the distribution 
of activity on a comet surface. A basic water sublimation model driven by solar inso- 
lation is used to derive total production rates for different nuclei shapes and spin axis 
orientations using the orbital parameters of 67P/Churyumov-Gerasimenko. The slopes of 
production rates versus heliocentric distance are calculated for the different model setups. 
The standard (homogeneous) outgassing model confirms the well known result regarding 
the heliocentric dependence of water production rate which remains invariant for differ- 
ent nuclei shapes as long as the rotation axis is perpendicular to the orbital plane. When 
the rotation axis is not perpendicular, the nucleus shape becomes a critically important 
factor in determining the water production curves as the illuminated cross-section of the 
nucleus changes with heliocentric distance. Shape and obliquity can produce changes in 
the illuminated cross section of up to 50% over an orbit. In addition, different spin axis 
orientations for a given shape can dramatically alter the pre- and post-perihelion produc- 
tion curves, as do assumptions about the activity distribution on the surface. If, however, 
the illuminated cross-section of the nucleus is invariant then the dependence on the above 
parameters is weak, as demonstrated here with 67P/Churyumov-Gerasimenko shape. The 
comets Hartley 2 and Wild 2 are shown to yield significantly different production curve 
shapes for the same orbit/orientation as 67P/CG, varying by as much as a factor of 3 as a 
result of only changing the nucleus shape. Finally, we show that varying just three basic 
parameters, shape, spin axis orientation, and active spots distribution on the surface, can 
lead to arbitrary deviations from the expected inverse square law dependence of water 
production rates near 1 AU. With the obtained results we cannot avoid the conclusion 
that, without prior knowledge of basic parameters (shape, spin-axis orientation, activity 
locations), it is difficult to reveal the nature of cometary outgassing from the heliocentric 
water production rates. Similarly, the inter-comparison of water production curves of two 
such comets may not be meaningful. 
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4.2 Introduction 


Comets are icy bodies which heat up as they approach the sun and sublimate volatile 
material from their nuclei to create a surrounding tenuous gas coma. Water vapor is 
generally the dominant component of the coma, especially around perihelion (Bockelée- 
Morvan, D. et al. 2002, Hässig et al. 2015), however at larger heliocentric distances where 
it is colder and the water ice is more stable, more volatile species such as CO can dom- 
inate the coma, as is the case for 29P/Schwassman-Wachmann (Bockelee-Morvan et al. 
2010). Biver et al. (1997) show that for comet Hale-Bopp, the outgassing transitions 
from CO to OH dominated at approximately 3.5 au and Laeuter et al. (2018) show that 
the water production rate can be lower than the production rates for CO and CO, from 
67P/Churyumov-Gerasimenko beyond 3.5 au post-perhelion. The water production rate, 
Q [molec. s^!], increases as the distance between the nucleus and the sun decreases due 
to increasing solar flux, nearly with an inverse square law dependence, but this is helio- 
centric distance dependent: for volatile ices such as CO and CO,, the production rates can 
be theoretically approximated by an inverse square law (r;?) at distances less than —4 au; 
for water, which has a higher sublimation temperature, the slope is usually much steeper 
between 1-4 au, and only approaches r;? within 1 au of the sun (Cowan and A’Hearn 
1979). As noted in Steckloff et al. (2015) and Steckloff and Jacobson (2016), by fully 
considering the temperature dependence of the sublimation coefficient (Gundlach et al. 
2011), the production rate approaches r;? within 0.5 au. 

Cometary water production rate curves are unique for each object and show few 
common features. Some comets exhibit reasonably symmetric water production rates 
about perihelion (Hale-Bopp, Kührt (1999)), while others show asymmetric behaviour 
(67P/Churyumov-Gerasimenko, Hansen et al. (2016)). In addition, the peak produc- 
tion does not necessarily occur at perihelion as it approximately does for 2P/Encke and 
46P/Wirtanen, but it can be offset to appear post-perihelion (as for 30P/Reinmuth 1 and 
6P/d' Arrest) or even pre-perihelion (as for 22P/Kopff and 81P/Wild 2) (Bondarenko and 
Medvedev 2011). 

The observed slopes for the water production rate also vary between comets. From the 
Rosetta mission, several attempts have been made to estimate the slopes for the change 
in the water production rate with heliocentric distance (Q — r") for 67P/CG (Hansen 
et al. 2016, Wedlund et al. 2016, Marshall et al. 2017), giving values of N = -5.3, -3.8 
and -7.1 pre-perihelion, and N = -7.1 and -4.3 post-perihelion. Values for N have also 
been obtained for other comets, including Hale-Bopp (N = -1.88, Russo et al. (2000)) and 
153P/Ikeya-Zhang (N = -3.21, Russo et al. (2004)). It is worth pointing out that all of 
these values have been obtained over different ranges of heliocentric distance. 

Although suitable from an observational point of view, we raise a question: what is 
the exact physical relevance of these slopes to the activity and given nucleus? And how 
can slopes for different comets be compared? 

In particular, we explore how the three basic factors - nucleus shape, obliquity of 
the rotation axis, and activity distribution - affect the water production rate of comets 
and their respective slopes!. It has been known for a long time (Sekanina 1981) that the 
orientation of the spin axis is very important for the light curves (brightness) of comets, 





'The Python code used in this chapter is available to download from: 
https://gitlab.com/david_marshall_mps/comet-outgassing-QvRh 
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with highly oblique objects like 6P/d' Arrest being able to sustain their brightness weeks 
after perihelion. We also discuss how our simplified model can be improved to add further 
realism. However, with more free parameters added to the problem it is unlikely that 
we can break the demonstrated degeneracy, if at least the basic nucleus shape, spin-axis 
orientation, and the active spot distributions are unknown. 


4.3 Sublimation model 


Following early works such as Watson et al. (1962), Cowan and A'Hearn (1979) and 
Weissman and Kieffer (1981), as well as Steckloff et al. (2015), the sublimation of water 
per unit area from an icy body can be written as: 


Sl = A,) S 


2 
T, 


€0T* + Z(T)Lices (4.1) 


where S, is the solar flux with a value of 1367 Wm’, A, is the bolometric Bond albedo 
with a value of 0.01 (Lamy et al. 2007), 7, is the heliocentric distance (dimensionless), € 
is the emissivity with a value of 0.9, o is the Stefan-Boltzmann constant, T is the surface 
ice temperature, Z(T) is the sublimation rate [kg s! m ?] and L;, is the latent heat of 
water ice with a value of 2.6x10° J/kg, taken to be a constant. This is similar to the form 
given in Keller et al. (2015). The rate of sublimation can be calculated from the equations 
(Panale and Salvail 1984, Langmuir 1913): 


Z(T) = 2P(T)/(avin(T)) (4.2) 
P(T) 2 3.56 x 10? exp(-6141.667/T) [kgm ' s ?] (4.3) 


v„(T) = vSRT /ny. (4.4) 


The constants in Eq. 4.3 are from solving the Clausius-Clapeyron equations. v; is the 
thermal velocity, R is the gas constant and u is the molar mass of water. 

The water production rate is simply the product of the sublimation rate of a region, 
Zregion(/n), and the sublimating area of the region, A,cgion. As we only use digital shape 
models in this work, the regions we describe are the triangular facets of a shape model. 
Summing over all facets, where N is the number of facets and i is the ith facet, gives the 
total water production rate: 


Olr) 7 X? Z(ri)A;. (4.5) 


For the sake of simplicity and to isolate the effects of spin axis orientation, shape, and 
the distribution of activity on the surface, Eq. 4.1 is presented only for highly idealized 
conditions of pure surface ice with zero thermal inertia. For the same reason we also 
neglect the temperature dependence of the sublimation coefficient (Gundlach et al. 2011). 
We assume that the cometary objects in this work are only made of water ice. These 
idealistic simplifications limit the validity of our model to heliocentric distances up to ~3- 
4 au (Meech and Svoren 2004). At larger distances, the heat conductivity term becomes 
significant and is dominated by the radiative part (Skorov et al. 2017). However, for 
detailed modeling, the composition, structure, and microphysical properties of the near 
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surface layer must be known. We acknowledge that cometary activity is more complex 
than assumed here (the so called Whipple model (Whipple 1950)), and additional process 
such as ice phase changes near the sub-surface, extended sources (such as inferred for 
103P/Hartley2), and super-volatile release impact water production rates and so does their 
surface distribution. Nevertheless, all of these remain poorly constrained apart from a 
few comets visited by spacecrafts (usually over a limited r, range) and would only add 
unnecessary complexity and additional unconstrained parameters. 

With the stated assumptions, when the total flux term (the left side of Eq. 4.1) is 
low, most of the energy goes into the thermal re-radiation term, and the sublimation term 
is negligible. As the total flux increases, more energy goes into sublimation and from 
approximately 50 - 1000 Wm”? (or 160 - 200 K surface temperature), the sublimation 
and thermal re-radiation terms on the right side are both similar in magnitude. At fluxes 
greater than 1000 Wm ? (or 200 K surface temperature), the sublimation flux accounts for 
more than 90% of the total flux and dominates over the thermal re-radiation term. When 
the flux is large (>1000 Wm ?), the thermal term does not contribute much, and Eq. 4.1 
reduces to: 


M x E Lo (4.6) 
replacing Z with Eq. 4.5. It can be seen that for large fluxes, for instance when a comet 
is close to the Sun, the water production rate follows an inverse square law for the helio- 
centric distance. 

For the numerical modelling described in this work, we use the sublimation model in 
Equation 4.1 to calculate the received flux and hence water production rate for a comet 
on an orbit around the Sun. We use a 67P-like orbit (eccentricity 2 0.64, semimajor axis 
— 3.46) and change the shape of the object on the orbit. For every facet of an object, the 
received flux can simply be calculated from the heliocentric distance and the orientation of 
the facet normal vector and the sun vector. Although important for very precise modeling 
of activity, shadowing and self-heating are not included in our considerations of real- 
shape illumination conditions. The received flux then corresponds to a temperature and 
sublimation rate given in a pre-calculated lookup table. The production rate for each 
facet is found by multiplying the facet sublimation rate by the facet area. The orbit is 
divided into fifty intervals and at each interval, the temperature, sublimation rate and 
production rate is determined for every facet over a single comet rotation. The average 
sublimation and production rates are then calculated for each heliocentric distance by 
averaging over all facets over one rotation. For simplicity we neglect the effects of self- 
heating for irregular shapes, which does not invalidate the main points of this work as 
demonstrated in Keller et al. (2015). 

The model does not need to account for a rotation rate. At each heliocentric distance 
interval, the comet is simply rotated and the water production is calculated for each inter- 
val during a comet day. 

This enables us to probe the effects of different shapes on the change in production rate 
with heliocentric distance. We also investigate the effect of the obliquity angle, defined as 
the angle between the rotation axis and the normal to the orbital plane, with the rotation 
axis in the same direction as the semimajor axis. In addition, we consider ®, the argument 
of the subsolar meridian at perihelion (hereafter, called argument ®) defined as the angle 
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swept out by the rotation axis from the semi major axis (Sekanina 1981). Obliquity and 
argument ® are rotations of the spin axis in two different perpendicular planes. We choose 
this method to demonstrate the effects of these two angles of axis orientation separately. 
Finally we study the effect of random and non-random distributions of activity on the 
surface. 


4.4 Results 


4.4.1 Effect of comet shape and obliquity 


We use six comet shapes (a sphere, 67P (Preusker et al. 2015), 81P/Wild2 (Farnham 
et al. 2005), 103P/Hartley 2 (Farnham and Thomas 20132), 1P/Halley (Stooke 2002) and 
9P/Tempel 1 (Farnham and Thomas 2013b)) and solve Eq. 4.1 for each comet around a 
67P-like orbit. This is done for obliquity values of 0°, 30°, 45°, 60° and 90°, keeping ® 
constant (equal to 0). 

We note that these shape models come from a variety of space mission flybys which do 
not have full coverage and only tens of metre resolution (except for the 67P shape model). 
However, this is not a problem for us as we only interested in how the shape affects the 
production rate curve. In addition, the errors and resolution of the shape models are not 
important for our model, the overall nucleus shape of the nuclei is the most important 
thing and will determine how the water production rate curve changes. 

The change in the normalized cross-section of the illuminated active area for each 
comet at each obliquity is shown in Fig. 4.1. Each curve is normalized to its maximum 
value. We define the cross-section of the active illuminated area, A, as a summation over 
all facets: 

A-LAy;cos(0y) (4.7) 


where A; is the area of a single facet and 6; is the angle between the facet normal vector 
and the sun-facet vector. If the angle is greater than 90°, it is set to 90° so that the flux is 
zero for non-illuminated facets. 

For the spherical comet, the illuminated area does not change with heliocentric dis- 
tance for any obliquity. In fact, regardless of shape, at 0° obliquity the cross-section of 
the active illuminated area is constant for any object with heliocentric distance. Despite 
67P’s irregular shape, its illuminated area does not change much, dropping by only 5% 
for an obliquity of 90°. 

When the rotation axis is at an oblique angle, the nucleus shapes of 81P/Wild 2 and 
9P/Tempel 1 cause significant changes in the illuminated area during an orbit. Addi- 
tionally, the larger the obliqueness of the rotation axis is, the larger the change in area. 
For both comets, the illuminated area decreases on approach to the sun before suddenly 
increasing again within | au from perihelion. 

Similarly, for 103P/Hartley 2 and 1P/Halley, the elongated shapes of the nucleus result 
in changes in the illuminated area, this time increasing on perihelion approach and then 
rapidly decreasing before closest approach. For these two comets, the (illuminated) area 
changes by more the 50% from maximum to minimum. 

The change in the illuminated area with heliocentric distance has a profound effect 
on the production rate curves, as seen in Fig. 4.2. The production rate for a sphere at 0° 
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obliquity is shown in dark blue, and the gradient of the production rate curve is shown in 
the bottom panels. A comet of any shape at 0° obliquity exhibits a very similar normalised 
production rate and slope, and a sphere at any obliquity also gives the same normalised 
result. 

The results for 67P at a rotation axis angle of 60° is in light blue. Even though we 
consider a large oblique angle, the production rate curve and the corresponding slopes 
do not change much in comparison to the spherical case with no axial tilt. Despite the 
shape and obliquity of 67P, it exhibits similar behaviour to a sphere, as a results of its 
coincidental invariance of illuminated cross-section with inclination. 

The dark and light orange lines lines in Fig. 4.2 represent the production rate and 
slopes for 81P/Wild 2 and 103P/Hartley 2 respectively, at 60° obliquity. In contrast to 
67P, these comets yield notable changes in production rates and gradients when obliquity 
changes, which is due to their shapes that do result in changes of total illuminated area. 

In the example for 103P/Hartley 2, the area starts to increase from large heliocentric 
distances until it reaches a maximum and then decreases in the final approach to perihe- 
lion. As a result, the production rate curve is steeper than the other cases far from the sun 
but less steep close to the sun, and the slope even becomes positive at perihelion. The 
increasing flux as the comet approaches the sun should increase the production rate but 
the decreasing illuminated area means that there is less flux available for sublimation. 

For Wild 2, the opposite behavior occurs. The illuminated area starts to decrease 
from large heliocentric distances until it reaches a minimum and suddenly increases as 
it approaches perihelion. The sudden increases in area results in steeper production rate 
slopes around perihelion than for any of the tested objects and obliquity setup. In this 
case, not only does the flux increase for the comet, but also the area, so the increase in 
water production is enhanced and the slope is steeper. 

It is therefore evident that, before even considering any properties of the comet nu- 
cleus, the water production rates and their corresponding slopes are heavily dependent 
upon the obliquity of the rotation axis and shape. 


4.4.2 Effect of activity distributions 


In the previous section, we considered all the surface of the comet to be active. This is 
unlikely to apply in reality. For example, Keller et al. (2015) and Marschall et al. (2016) 
show that only certain local activity can explain the outgassing observations of 67P. In this 
section, we show how features in the production rate curves and their slopes are affected 
by random and non-random distributions of active regions on the comet. 

We examined the production rate curves for six different distributions of activity on a 
sphere: 100% active, 50% randomly distributed activity, 20% randomly distributed activ- 
ity, an active region on the northern hemisphere, an active region on the southern hemi- 
sphere, and an active belt at 0° longitude. 

For an obliquity of 0°, all the distributions produced similar production rate curves; 
only the absolute value of Q is affected. The illuminated area stays constant with helio- 
centric distance and only the magnitude changes depending on how much of the surface 
is active. 

For an obliquity of 60° though, there are significant differences between the produc- 
tion rate curves for each kind of distribution. These are shown in Fig. 4.3. In the two 
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Figure 4.1: The change in the normalised cross-section of the illuminated area for six 
selected comet shapes (sphere, 67P, 81P, 103P, 1P and 9P) at five different obliquities (0°, 
30°, 45°, 60° and 90°) with distance in au from perihelion (on a 67P-like orbit, perihelion 
occurs at 1.24 au). 
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Figure 4.2: Normalised water production rate and slopes of different shaped comets 
(sphere, 67P, 81P, 103P). A sphere with an obliquity of 0° is shown as a reference. Any 
shape at 0° has a similar production rate curve. Comets 67P, 81P and 103P are shown 
at obliquities of 60°. The left and right panels correspond to pre- and post-perihelion re- 
spectively. Pre-perihelion distances are denoted as negative and post-perihelion distances 
are given as positive. For a 67P-like orbit, perihelion occurs at 1.24 au (or -1.24 au in the 
pre-perihelion figures on the left). 
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Figure 4.3: Normalised area, water production rate and slopes of spherically shaped 
comets with different distributions of activity (100% active, 50% active, 20% active, ac- 
tive region on the north hemisphere, active region on the southern hemisphere, active belt 
at 0° longitude) at 60° obliquity. 


randomly distributed cases, the illuminated areas are reasonably constant although they 
can have fluctuations based on how the random activity is distributed. As expected, these 
fluctuations do not change the production rate curves much and slopes are comparable to 
the fully active case. 


When considering an active region on the northern side, the illuminated area increases 
dramatically as the comet approaches perihelion, but remains on the night side for long 
period going to aphelion, due to the assumed obliquity of the rotation axis. The production 
rate slope starts out very steep, but around perihelion, the curve and slope are comparable 
to the randomly distributed cases. When the active region is on the southern side, the 
opposite happens. The active region faces the sun at aphelion and is away from the sun at 
perihelion. Hence, the illuminated area decreases with heliocentric distance and the slope 
is therefore less steep at perihelion and becomes positive. The peak in production is offset 
before and after perihelion. 


In the final case, activity is confined to a belt at 0? longitude. Here, the illuminated 
area decreases until it reaches a minimum and then increases towards perihelion. This is 
identical in behaviour to a fully active comet shaped like Wild 2, and both of their produc- 
tion rate curves behave in the same way, with a significant steepening at perihelion. This 
shows that a sphere with a defined region of activity can have the same production rate 
curve as an irregularly shaped fully active comet. As a result, it is difficult to disentangle 
the effects of shape and activity from only the production rate curves. 
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4.4.3 Effects of obliquity and ® 


In all the simulations so far, the pre- and post-perihelion curves have been symmetrical 
around perihelion. This is because the obliquity angle has always been orientated towards 
perihelion (parallel to the semi-major axis), so the change in illuminated area during ap- 
proach and recession is symmetrical. We now introduce the argument ®, as defined in 
Section 4.3. Its effects are actually well known in dealing with light curves of asteroids 
and comet brightness observations. Here we show that variations in also have a strong 
effect on the production rate curves. 


We show four selected examples in Fig. 4.4 to illustrate the main point for a 103P/Hartley 
2 shaped comet: an obliquity of 0? and an argument ® of 60°; an obliquity of 60° and an 
argument ® of 90°; an obliquity of 60° and an argument ® of 60°; and an obliquity of 80° 
and an argument ® of 10°. 


For the first case, there is no change in the illuminated area. ® changes the angle 
between the rotation axis and the semimajor axis, but as the obliquity angle is zero, the 
rotation axis is still perpendicular to the orbital plane. As a result, the illuminated cross 
section does not change with heliocentric distance and the production curves remain un- 
altered. 


With an obliquity of 60° and argument ® of 90°, the change in illuminated area is 
the inverse of the change for 103P/Hartley 2 in Section 4.4.1, which had an obliquity 
of 60°, too. Argument ® orients the rotation axis away from the comet-sun vector at 
perihelion/aphelion, rather than in the same direction, as in Section 4.4.1. Now, the area 
decreases with heliocentric distance before increasing at perihelion. This has the effect of 
creating a steep production rate curve around perihelion, similar to the case of Wild 2 in 
Fig. 4.2. Importantly, differently shaped comets can be seen to have similar behaviours in 
their production rate slopes for different values of ®. The pre- and post-perihelion curves 
are still symmetric around perihelion though. 


Asymmetry in the water production can be produced for other values of obliquity and 
argument ®. With the obliquity at 60° and the argument ® at 60°, a slight offset can be 
seen, moving the peak production from 1.24 au (perihelion) to 1.26 au post-perihelion, 
equivalent to a delay of approximately 16 days. The increase in area post-perihelion is 
large enough to compensate for the loss of flux from moving out to larger heliocentric 
distances. 


Finally, for case of an obliquity of 80° and ® of 10°, peaks in the production rate can 
be seen pre- and post-perihelion as the illuminated area reaches a maximum before and 
afterwards. The peak production can now be seen at 1.4 au pre- and post-perihelion, shifts 
of about 50 days. 


Our modelling is incapable of producing an offset in water production for 67P for 
an obliquity of 52° and argument ® of 20°, the estimated values for 67P. This is due to 
the nearly constant behaviour of the illuminated cross section with heliocentric distance. 
There is not a sufficient increase in area to overcome the decreasing amount of flux re- 
ceived as heliocentric distance increases. A defined active region would have to be used 
to create an offset in our model. 
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Figure 4.4: Normalised area, water production rate and slopes of 103P/Hartley 2 shaped 
comets with different obliquities and ® angles (obliquity = 0°, ® = 60°; obliquity = 60°, 
® = 60°; obliquity = 80°, ® = 10°). The pre- and post-perihelion curves are plotted on 
the left and right panels respectively. 


4.5 Conclusions 


We applied a simplified model of water outgassing rate for different cometary bodies to 
investigate how cometary nucleus shape, spin axis orientation and activity distribution 
determine the observable dependence with heliocentric distance. We acknowledge that 
our modeling relies on strong simplifications and does not correspond to any realistic 
cometary object. Nevertheless, even in this idealistic form it is sufficient to demonstrate 
the importance of the investigated parameters on the heliocentric production rates of water 
within 3 au (ie, spin-axis orientation, shape effects, and activity distribution). The appli- 
cation of more comprehensive and detailed models (accounting for self-heating, shadow- 
ing, gas diffusion through porous dust mantles with moving boundary, depth variation 
of conductivity, micro-fracturing, re-condensation effects, perhaps ice-crystallization and 
number of other effects) is beyond the scope of this work, but they clearly merit future 
investigations for objects where at least the basic nucleus shape and spin-axis orientation 
are known. In the introduction we asked: Is there a qualitative and/or quantitative way 
to interpret the behaviour of Q(r;) of a comet (for which the exact shape, spin axis and 
activity distribution are not known), and how to extract meaningful inter-comparisons of 
such observations for different comets? 
The following general conclusions follow from the presented results: 


e For homogeneous, solar driven activity, the illuminated cross-section of the nucleus 
as it travels on its orbit governs the water production rate, Q(r;), and its gradient 
with heliocentric distance. 


e This variability, even under an assumption of homogeneous icy surface, is signifi- 
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cantly modulated by the basic parameters: nucleus shape and spin axis orientation. 
For example, a comet shaped like 103P/Hartley 2 with a 90° obliquity produces 
changes in the illuminated cross-section of 50% compared to 0° over an orbit, 
whereas changes of less than 20% are seen for a comet shaped like 9P/Tempel 1 at 
90° obliquity. 


e The often assumed (approximately) accurate Q(r;) ~ r ? is valid only for homo- 
geneously active objects at zero obliquity. This relationship holds only over a very 
narrow heliocentric range around 1 au. 


e [tis insufficient to derive a single value for the slope in the form, Q(r;) = Ary B over 
a large range of heliocentric distances (covering data beyond 3 au), as the slope, B, 
is also a function of distance. Even in the simplest case of a sphere with 0 degree 
obliquity, the production rate slope varies from a value of -4 at 3 au to -2.5 at 2 au. 
Therefore, inter-comparisons of such values among different comets may not be 
informative. 


e [f there is a non-uniform distribution of active regions on a comet body (possibly 
due to any combination of dust cover/mantle, extended sources, albedo, etc), vir- 
tually any shape of O(r,) can be obtained. This is especially true for irregularly 
shaped objects with non-zero axis obliquity and ® angle. 


e The axis orientation (obliquity and ®) is a sufficient condition for irregular shapes 
to yield asymmetric Q(r;) and its slopes. The strength of this effect is modulated 
by the nucleus shape (illuminated cross-section variations). 


Inferring physical properties of (unresolved) comets from their production rate curves 
is an exceptionally difficult task and caution should be exercised when interpreting these 
curves. 

The above conclusions have strong implications on inter-comparisons of the Q(r;) for 
different comets, especially when nucleus shape and spin axis orientation are not known. 
We argue that direct comparison of slopes of heliocentric production rates obtained from 
different distances as well as mixing pre and post-perihelion data make it very difficult 
to draw consistent conclusions. It is not rare that such comparisons are published in 
the literature (usually a result of observational constraints), with one recent example in 
Wedlund et al. (2016). Such comparisons, as shown in this work, cannot be physically 
meaningful except by chance. 

At last, a suitable example where the shape and spin-axis orientation is well deter- 
mined but the slope of Q(r;)[H5O] do not obey the r;? approximation is the comet 67P. 
This is surprising, since it retains nearly a constant cross-section despite its spin axis 
orientation throughout its orbit as shown. From Rosetta measurements the production 
rate slope ranges from -3 to -8 (Hansen et al. 2016, Marshall et al. 2017), which based 
on our results can be explained only by an inhomogeneous distribution of water activity. 
Observations from MIRO and ROSINA have already implied that the activity is not ho- 
mogeneously distributed, with activity in particular regions being stronger than in other 
areas (Marshall et al. 2017, Marschall et al. 2016, Fougere et al. 2016). The observed de- 
lay in peak production rate for 67P (Hansen et al. 2016, Marshall et al. 2017, Bondarenko 
and Medvedev 2011) must be a combination of the factors described here, or additional 
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processes not considered. Applying more sophisticated models to detailed observations 
of 67P may provide a better fit, but to what degree such fit is unique is still open. This 
fact only underlines our general argument using a simple model. 
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This chapter has been published in A&A as the journal article: Marshall, D., et al. "Ther- 
mal inertia and roughness of the nucleus of comet 67P/Churyumov-Gerasimenko from 
MIRO and VIRTIS observations." Astronomy & Astrophysics 616 A122 2018. Repro- 
duced with permsission © ESO 


5.1 Summary 


Using data from the Rosetta mission to comet 67P/Churyumov-Gerasimenko, we evaluate 
the physical properties of the surface and subsurface of the nucleus and derive estimates 
for the thermal inertia (TI) and roughness in several regions on the largest lobe of the 
nucleus. 

We have developed a thermal model to compute the temperature on the surface and in 
the uppermost subsurface layers of the nucleus. The model takes heat conduction, self- 
heating, and shadowing effects into account. To reproduce the brightness temperatures 
measured by the MIRO instrument, the thermal model is coupled to a radiative transfer 
model to derive the TI. To reproduce the spatially resolved infrared measurements of the 
VIRTIS instrument, the thermal model is coupled to a radiance model to derive the TI and 
surface roughness. These methods are applied to Rosetta data from September 2014. 

The resulting TI values from both instruments are broadly consistent with each other. 
From the millimetre channel on MIRO, we determine the TI in the subsurface to be 
«80 JK^!m^?s^9? for the Seth, Ash, and Aten regions. The submillimetre channel im- 
plies similar results but also suggests that higher values could be possible. A low TI 
is consistent with other MIRO measurements and in situ data from the MUPUS instru- 
ment at the final landing site of Philae. The VIRTIS results give a best-fitting value of 
80 JK^! m ?s ?? and values in the range 40-160 JK^! m?s ?? in the same areas. These 
observations also allow the subpixel scale surface roughness to be estimated and com- 
pared to images from the OSIRIS camera. The VIRTIS data imply that there is significant 
roughness on the infrared scale below the resolution of the available shape model and 
that, counter-intuitively, visually smooth terrain (centimetre scale) can be rough at small 
(micrometre - millmetre) scales, and visually rough terrain can be smooth at small scales. 
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5.2 Introduction 


Comets are considered to be remnants of the pristine material which formed during the 
early solar system, and as such, they may hold clues to their local formation environment 
and the evolution of the solar system. One way to assess the pristine nature of a comet 
is to measure the thermal properties of the surface layer because these control the energy 
transport through the surface and thus how much material is still unprocessed by solar 
insolation since its formation. 

Energy transport through the surface of the nucleus is a very complex process (e.g. 
Gortsas et al. 2011). Itis mainly initiated by the incoming solar insolation and the dust-ice 
matrix of the surface allows for the propagation of energy into the subsurface layers. The 
dust and ice can emit radiation and the volatile ice can also sublimate from the cometary 
material when heated. Since the volatiles are thought to be distributed throughout the 
surface and subsurface layers (De Sanctis et al. 2015), sublimation of ices could arise from 
different depths within the nucleus as heat is transported internally. Dust can then also be 
expelled from the surface as a result of ice sublimation. In addition, the inhomogeneous 
surface topography leads to shadowing and self heating (Keller et al. 2015). 

Thermal inertia (hereafter TI) and surface roughness are two of the key mediators for 
the propagation of energy into a comet nucleus. The TI describes how efficiently the solar 
energy propagates into the surface and represents the ability of the surface to respond 
to the temperature forcing provided by the solar illumination. Roughness can locally 
enhance or reduce the amount of energy received and emitted at a given location. These 
two properties can be determined from the spatial and temporal surface temperature. 

The Rosetta mission was launched in 2004 and arrived at its target destination, the 
comet 67P/Churyumov-Gerasimenko (hereafter 67P), in August 2014. The study of 67P 
represents an opportunity to better understand the physical properties of the nucleus as 
inferred from the energy transport in a cometary nucleus. For two years, the spacecraft 
observed the behaviour of the comet as it approached and receded from the sun. The 
mission ended in September 2016 when the spacecraft descended onto the surface. There 
were 11 instruments onboard Rosetta, and we used data obtained by two remote sensing 
units to determine the thermal properties of the surface and subsurface: MIRO (Gulkis 
et al. 2007, Microwave Instrument for the Rosetta Orbiter) and VIRTIS (Coradini et al. 
2007, Visible and InfraRed Thermal Imaging Spectrometer). These are described in the 
next section. 

Before the Rosetta mission, the Deep Impact spacecraft measured the thermal proper- 
ties of comets 9P/Tempel 1 and 103P/Hartley 2 (Groussin et al. 2013). The TIs for these 
two comets were found to be lower than 250 JK^!m^?s ?? for Hartley 2 and lower than 
45 JK !m?s 9? for Tempel 1 (30 upper limit). The regions with exposed water ice had 
temperatures over 100 K higher than the sublimation temperature (~200 K) , indicating 
the presence of ice-free dust at the subpixel scale. In addition, work by Davidsson et al. 
(2015) demonstrated that disc-integrated thermal emission from comets depends upon the 
type of roughness on the surface as well as the TI. 

Regarding 67P, in a previous work by Schloerb et al. (2015), the TI in the Imhotep and 
Ash regions was calculated from MIRO data. These authors derived a low TI between 10- 
30 JK^! m^?s ?? and found that the MIRO millimetre (hereafter mm) emission arises from 
a depth of approximately 4 cm, whilst the submillimetre (hereafter smm) emission comes 
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from approximately 1 cm. Their work also shows residuals that are correlated to surface 
features, suggesting that regional heterogeneity may affect the surface properties of the 
comet. Choukroun et al. (2015) also calculated TI from MIRO data and obtained values 
of 10-40 JK^! m^?s?? and 20-60 JK^! m ?s 9? for the submillimetre and millimetre chan- 
nels, respectively. In addition, Gulkis et al. (2015) found low TIs of 10-50 JK^! m ?s 9? 
from MIRO observations as Rosetta first approached 67P. 

From the Philae lander, MUPUS (Spohn et al. 2007, Multipurpose Sensors for Sur- 
face and Sub-Surface Science) measured a TI of 85 € 35 JK^!m ?s ?? at the final landing 
site, Abydos (Spohn et al. 2015). This represents the best-fitting value from their thermal 
models to the infrared radiometer data from the MUPUS thermal mapper. This is higher 
than the MIRO measurements but this could be attributed to heterogeneities in the surface 
layer, such as porosity, temperature, dust-to-ice ratio, composition, dust distribution, and 
layer thickness (Spohn et al. 2015). The MUPUS measurement is also larger than unre- 
solved observations made by the Spitzer infrared space observatory, which gives a low TI, 
20 JK!m^?s 9? (Lamy et al. 2008), and is more consistent with the results obtained from 
MIRO measurements. The results from thermophysical modelling applied to light curves 
also gave a low TI value equal to 15 JK^!m^?s ?? (Lowry et al. 2012). 

For the asteroid Lutetia, which Rosetta passed during its flight to 67P, Keihm et al. 
(2012) used observations obtained by MIRO and VIRTIS to determine the TI and rough- 
ness. These authors found that a low TI («30 JK ! m ?s 95), which increased with depth, 
and a 5046 fractional surface coverage with hemispherical craters was required to fit the 
model to the observations. 

We aim to determine thermal properties of the surface layer of the nucleus of comet 
67P, such as the TI and roughness, by analysing data from the VIRTIS and MIRO experi- 
ments. Both of these datasets provide measurements of the emitted and reflected radiance 
of the cometary nucleus. In the case of VIRTIS, the emitted radiance is related to the sur- 
face temperature (uppermost few tens of microns), whereas for MIRO, which operates at 
longer wavelengths, the emitted radiance is related to the subsurface temperature profile. 
We focussed our analyses on data obtained in September 2014. This choice was motivated 
by the opportunity to compare our final results with previous publications from the MIRO 
science team, such as by Schloerb et al. (2015), Gulkis et al. (2015), and Choukroun et al. 
(2015), who analysed MIRO data obtained during the same time period although from 
different regions. 

In the next section, we describe the MIRO and VIRTIS instruments. Then, the first 
step towards estimating the TI involves assessing the datasets to find suitable data for 
comparison (Section 5.4.1). The thermal model applied to the data from both instruments 
and used to calculate the temperature gradients is described in Section 5.4.2, and the 
details of the MIRO radiative transfer model are given in Section 5.4.3. The VIRTIS 
modelling is described in Section 5.4.4. The results from MIRO and VIRTIS are given in 
Sections 5.5.1 and 5.5.2, respectively. 


5.3 Instruments 


The MIRO instrument was a microwave spectrometer consisting of a 30 cm offset parabolic 
reflector telescope; a millimetre heterodyne receiver, operating at a centre-band frequency 
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of 188 GHz (1.6 mm wavelength); and a submillimetre heterodyne receiver, operating at 
562 GHz (0.5 mm wavelength) (Gulkis et al. 2007). The frequency bands of the millime- 
tre and submillimetre receivers had angular resolutions of 23.8 arcmin (6.9 mrad) and 7.5 
arcmin (2.2 mrad), respectively, and each contained a single broadband channel for the 
measurement of near-surface temperatures. Additionally, a Chirp Transform Spectrome- 
ter (Hartogh and Hartmann 1990, CTS) was connected to the submillimetre receiver for 
the observation of absorption and emission lines of water, carbon monoxide, methanol, 
and ammonia. The work presented here focusses on the temperatures measured by the 
two broadband channels. 


The instrument observed the source for 30 seconds and was calibrated in situ every 
34 minutes by observing onboard hot and cold calibration targets. In September 2014, 
these typically had temperatures of 18°C and -47°C. The output from the receivers was 
proportional to the observed power and converted to an antenna temperature, T4, by 


TA = P/k, (5.1) 


where P is the observed power density and k is the Boltzmann constant. In this work, we 
converted the measured antenna temperatures from MIRO to a brightness temperature, Tg, 
which is defined as the required temperature of a black body, which fills the antenna beam 
pattern to produce the observed power. A Planck function was used to calculate Tz. We 
assumed that each beam could be approximated as a Gaussian beam, and a multiplicative 
correction factor of wa was applied to the brightness temperatures from the submillimetre 
channel as a result of the difference in beam efficiency between the two channels, as 
described by Schloerb et al. (2015). 


The VIRTIS instrument was a hyperspectral imager that consisted of a high-spectral- 
resolution point spectrometer (VIRTIS-H, 1.88—5.0 um) and two mapping channels 
(VIRTIS-M-VIS, 0.22-1.05 um; VIRTIS-M-IR, 0.95-5.1 um) (Coradini et al. 2007). We 
used data acquired by the IR mapping channel VIRTIS-M-IR (hereafter abbreviated as 
VIRTIS), the cryocooler of which failed in May 2015, thereby ending its science out- 
put. Its instantaneous field of view was 250 urad x 250 urad. At each basic acquisition, 
VIRTIS recorded a frame, a spectrally resolved spatial line (432 spectral bands times 256 
spatial samples; 9.5 nm spectral sampling). A series of consecutive frames forms a data 
cube, a spectrally resolved two-dimensional image. For the cubes considered here, a basic 
acquisition typically had an exposure duration in the order of seconds, whereas the acqui- 
sition of a cube takes of the order of tens of minutes up to hours, a time span over which 
the observation and illumination geometry can change considerably from comet rotation 
and spacecraft motion. 


The VIRTIS infrared spectrometer was mostly sensitive to the surface temperature 
(uppermost few tens of jum), rather than the centimetre subsurface as in the case of MIRO. 
The longest wavelengths that VIRTIS could measure generally provide the best diagnos- 
tics since the contribution from reflected flux is less important at these wavelengths. We 
always see the rising wing on the short wavelength side (the Wien’s part) of the thermal 
emission peak, starting from about 3.0-3.6 jum (Filacchione et al. 2016). 
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5.4 Methods 


5.4.1 Observational overlap 


In order to better determine the thermal properties of the nucleus in some regions, we 
decided to combine MIRO and VIRTIS data, providing more observational constraints 
than with a single instrument and enabling us to investigate possible heterogeneities be- 
tween surface and subsurface thermal properties. We therefore identified when MIRO and 
VIRTIS were observing the same location on the comet at the same time. The data were 
averaged into one-minute time intervals and using a digital shape model of the nucleus, 
we could determine the footprints on the nucleus of the two instruments at a given time. 
We used a decimated version of shape model "cg-dir. spg-shap7-v1.0", which has 124938 
facets (Preusker et al. 2017) with a horizontal resolution of approximately 20 m, and the 
SPICE software (Acton 1996). This way, we could determine the averaged viewing ge- 
ometry of each observation from the positions of the spacecraft, the comet, and the sun, 
and then calculate the viewing angle (angle between facet normal vector and the vector 
connecting the facet and the spacecraft) and beam size of the submillimetre and millime- 
tre beams, as well as the footprint of the observation and the facet illumination. Most of 
the time MIRO and VIRTIS observed different regions of the nucleus but occasionally 
they both looked at the same location. 

The two instruments acquired data in slightly different ways. The MIRO millime- 
tre and submillimetre beams probed a single circular area (assuming a Gaussian beam 
shape), giving two subsurface temperatures, one from each beam. At 30 km from the 
nucleus the millimetre and submillimetre beams had approximate diameters of 210 m and 
64 m, respectively. The millimetre beam hence observed between 55 and 367 facets of 
the decimated shape model for the observations given in Table 5.1, and between 3 and 
26 facets in the submillimetre beam. The VIRTIS instrument, on the other hand, was a 
slit spectrometer, observing a narrow swath across the surface at a single exposure. At 
30 km, it instantaneously observed an area of approximately 8 m x 2000 m and could 
resolve the temperature of the 256 facets placed along the swath. We therefore looked 
for observations in which the MIRO beam intersects the VIRTIS band. Since both instru- 
ments operated at different wavelengths and spatial resolutions, the roughness is treated 
differently in each case. The MIRO instrument takes large-scale topographical roughness 
at the scale of the shape model (metres) into account when inferring the TI, and VIRTIS 
models subfacet roughness, since it has a higher spatial resolution. 

In September 2014, observations from the MIRO and VIRTIS instruments overlapped 
spatially and temporally ~150 times, and most of these occurred between September 1-2 
and September 12-15. For this work, we chose to focus on five of the best times; these 
times were selected because they all have well-constrained MIRO and VIRTIS observa- 
tions, good viewing geometries, and cover different areas of the nucleus. These times 
(specified in UT) are given in Table 5.1 along with the corresponding observed MIRO 
brightness temperatures. These observations provide data on five areas located on differ- 
ent geomorphological areas (El-Maarry et al. 2015). On September 2 and 15, 2014, the 
instrument probed areas located in the Aten region, while on September 1 and 13, it ob- 
served areas located in the Seth and Ash regions, respectively. On September 12, MIRO 
and VIRTIS probed an area located at the border between the Ash and Aten regions. The 
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Date Time smm mm Local solar S/C Beam Beam 
Ts Ts time distance latitude — longitude 
K K (hours) (km) (degrees) (degrees) 
Sep 1,2014 23:48 188 173 9.9 52 20 -143 
Sep 2, 2014 04:23 170 149 10.4 53 10 91 
Sep 12,2014 23:28 157 156 8.9 30 45 111 
Sep 13,2014 09:33 165 163 8.5 30 37 174 
Sep 15,2014 02:51 171 145 10.6 30 7 86 





Table 5.1: Brightness temperatures, Tg, from the submillimetre and millimetre channels 
for each MIRO observation, along with the local solar hour, spacecraft distance, latitude, 
and longitude at the centre of the beam footprint. With an assumed beam efficiency of 
0.94, we use errors on the brightness temperatures of 4 K. 


measurements correspond to the pre-landing mapping phase shortly after arrival at the 
comet. During this epoch Rosetta was orbiting 67P at a distance of 30 km, at 3.38 AU 
from the sun. 


5.4.2 Thermal model 


To interpret the MIRO and VIRTIS data, we need to calculate the temperature of the 
surface and subsurface in response to insolation. The thermophysical temperature pro- 
file of a comet nucleus can be determined by solving the heat conduction equation with 
the assumption that the heat induced by the solar radiation is transported into depth and 
emitted into space (Groussin et al. 2013). This can also be used to sublimate volatiles as 
considered by Capria et al. (2017), but we neglected this part in this work, since water 
sublimation is negligible at 3.38 AU. The subsurface temperatures can be approximated 
with a one-dimensional thermal model using the heat equation 


TZA &T(z,t) 

a m” 
where p is the bulk density in units of [kg m7], c, is the specific heat [J kg ! K^!], and 
K is thermal conductivity [W m~! K^!]. The partial derivatives areo and aren are for 
the changes of temperature with time, t, and depth, z, respectively. In Eq. (5.2), constant 
surface properties are assumed that do not depend on depth or temperature or time, and 
there is no internal heat source. The TI is defined as 








PCp 


(5.2) 


I -4 ipe, (53) 
where 7 has SI units of [J K! m? s ??]. This is the quantity of interest, and with these 
equations, the temperature is found as only a function of TI, assuming Bond albedo and 
emissivity are known. We take the density to be 532 kg m ? (Jorda et al. 2016) and assume 
the heat capacity to be 500 J kg ! K^!, estimated from Robie et al. (1970). 
For the surface boundary, the following condition is used to solve the heat equation 
S(1—A;, . OT 
SUA) cosi 2 eo T^ — k—— |.-o, (5.4) 


2 
Oz 
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where S, the solar flux at 1 AU, is equal to 1370 W m^, Ay is the Bond albedo, taken 
to be 0.0108, from a standard geometric value of 0.04 and phase integral value of 0.27 
(Lamy et al. 2007), e is the infrared emissivity, taken to be 0.95 (Birkebak 1972), o is the 
Stefan- Boltzmann constant, r, is the heliocentric distance, and i is the angle of incidence 
of the sun to the local surface. 

As a first step, the illumination geometry and projected shadows are calculated for 
each facet of the shape model over one complete nucleus rotation. This is carried out for 
three dates for which we have overlapping VIRTIS and MIRO data: September 2, 2014, 
September 13, 2014, and September 15, 2014. 

In the second step, we calculate the temperature of each facet with a thermal model 
that takes into account insolation, emission, self-heating, and heat conductivity with dif- 
ferent values of TI (see Eqs. (5.2 - 5.4) and Groussin et al. (2013)): 5, 10, 20, 40, 80, 160 
and 320 JK-!m~*s~°>. In this work, we test the TI at these discrete values. In the conclu- 
sions, rather than discuss error bars, we consider a range of acceptable values denoted by 
the tested values. The thermal model is run over several nucleus rotations for each date 
(between 8 and 22 rotations depending on the thermal inertia) with a time step of 12.4 
sec and up to a depth of 10 diurnal skin depths (between 2.2 cm to 1.4 m depending on 
the TI). Running to these depths ensures convergence of the temperature profiles up to 
ten wavelengths of the MIRO millimetre channel (1.6 cm), and by checking this conver- 
gence, we negated the need for a full seasonal model for the uppermost layers of interest. 
The temperature at the lowest depth was assumed to be 30 K, and a depth is used such 
that a = 0. To account for the large number of facets (124938 facets) and the physical 
processes such as heat conductivity and self-heating, the code has been optimized and 
parallelized. On a computer with four CPU cores, it takes approximately 48 hours to run 
a complete set of TI values at a given date. 

As an output, we obtain, for a given date, the temperature profile for each facet of the 
shape model from the surface to a depth of 10 MIRO mm wavelengths or more. Figure 5.1 
shows an example of temperature as a function of depth, for a TI of 80 JK~'m~’s~°°, on 
September 15, 2014 (UT 00:00). As expected for the dayside, the temperature decreases 
with depth. Figure 5.2 shows an example of surface temperature, for varying values of the 
TI in the range 5-320 JK ^! m ?s 9^, on September 15, 2014 (UT 00:00). The maximum 
temperature decreases with increasing TI (since more energy penetrates into the nucleus 
by heat conductivity), as does the diurnal thermal amplitude. 


5.4.3 Radiative transfer model for MIRO data 


The MIRO radiometer measures antenna temperatures at millimetre and submillimetre 
wavelengths, described in Section 5.3, from which brightness temperatures for the ob- 
served areas are derived. Using the thermal outputs from Section 5.4.2 as inputs to a ra- 
diative transfer model, simulated brightness temperatures (hereafter SBT) are compared 
to measured brightness temperatures (hereafter MBT) for the different MIRO obervations. 
The TI of the material layer is constrained if the SBT curve for a given TI intersects the 
MBT curve. The SBT varies as a function of the TI and the relative complex permittivity 
(hereafter permittivity) of the material taken into account in the radiative transfer. 

The permittivity, €, consists of a real part, €’, and an imaginary part, €", which describe 
the polarisability of a material and the capacity of motion of free charges in the material 
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230K 





115K 
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Figure 5.1: Temperature as a function of depth (on the surface and 1 cm, 2 cm, and 4 cm 
below the surface) for a TI of 80 JK^!m^?s^9?, for September 15, 2014 (UT 00:00). 


when an electric field is applied, respectively. Their values are relative to that of a vacuum, 
which is equal to 8.85x10^? F m*!. The real part of the permittivity is often referred to as 
the dielectric constant. The permittivity depends on the frequency, related to polarization 
mechanisms, composition, temperature, and bulk density of the material (Campbell and 
Ulrichs 1969, Mattei et al. 2014, Brouet et al. 2016). 

For a given observed area and at a given wavelength of observation, the SBT is deter- 
mined from the inverse Planck function taking into account the total spectral brightness 
intensity received from all facets covered by the MIRO beam. The spectral brightness 
intensity of a given facet varies as a function of the temperature profile of the subsurface 
and the permittivity of the material. It also depends on the emission angle, which refers 
to the orientation of the facet with respect to the main beam of the MIRO antenna. In this 
way, roughness at the scale of the shape model (of the order of metres) is treated in the in- 
terpretation of the MIRO measurements. A weighting function is applied to calculate the 
spectral brightness intensity of each facet. This function is based on an assumed Gaussian 
beam pattern for the MIRO instrument (Gulkis et al. 2007), where the centre of the beam 
contributes more strongly than the edges. 

For one facet, the spectral brightness temperature intensity BT, as a function of the 
effective spectral brightness intensity BT, and the surface emissivity, Y, is as follows 


BT) facet(@c) BT cf ¢(Oe) x YG, 6e), (5.5) 
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Figure 5.2: Temperature on the surface for various values of the thermal inertia in the 
range 5-320 JK^! ms 95, for September 15, 2014 (UT 00:00). 


where v is the frequency of the observation, 0,, the facet emission angle defined as the 
angle between the facet normal and the main beam of MIRO antenna during the observa- 
tion, and 6,, the transmission angle in the subsurface layers of the facet. The transmission 
across a specular boundary is considered in order to compute the emissivity of the facet 
surface using the Fresnel reflectivity coefficients. The effective spectral brightness in- 
tensity is computed integrating the contribution of each subsurface layer, weighted by a 
radiative transfer function, which expresses the extinction of the spectral brightness in- 
tensity as it propagates towards the surface. The effective spectral brightness intensity 
computed for a given facet, assuming a scatter-free homogeneous layer under conditions 
of local thermodynamic equilibrium, is expressed as follows (Ulaby et al. 2014) 


1 H — 
BT, ¢f@.) = -r X f B(v, T(z))e "C V1-sim 62/ex&qz — (5.6) 
V1 — sin(0,)/e x 6, 0 


where B(v, T (z)), the Planck function for a given frequency of observation and tem- 
perature, T(z), which varies with depth and is derived from the thermal model, and 6ó,;, 
the electrical skin depth of the layer at a given wavelength. The value H is the thickness 
of a layer equalling ten thermal skin depths. The electrical skin depth is inversely propor- 
tional to the imaginary part of the permittivity and proportional to the wavelength and the 


87 


5 Paper III - Thermal inertia and roughness 





square root of the dielectric constant. It characterizes how deep an electromagnetic field 
can penetrate into the material. 

Regarding the MIRO observations, we expect the submillimetre receiver to be sensi- 
tive to the thermal radiation integrated over the top subsurface of the layer and the mil- 
limetre receiver to be sensitive to the thermal radiation integrated over a greater depth. 
The electrical skin depth is constrained by a permittivity model based on experimental 
results obtained with lunar regolith samples reported in Olhoeft and Strangway (1975). 
The temperature dependence of the dielectric constant is considered to be linear (Calla 
and Rathore 2012, Brouet et al. 2015). 


5.4.4 Radiance model for VIRTIS data 


With the VIRTIS data, we focus our study on only two specific dates in Table 5.1: Septem- 
ber 12 and September 15, 2014. 

Fig. 5.3 shows the shape model of 67P rendered from the observing view of the space- 
craft at those two epochs, and the radiance measured at 4.95 um by VIRTIS around that 
time. This wavelength was chosen because it is the longest in the VIRTIS range that is 
not directly at the detector edge and is outside certain spectral features that may be real 
or else caused by calibration issues (Filacchione et al. 2016). The measured radiance is a 
combination of directional thermal emission and reflected solar radiation. The observed 
areas cover the regions of Aten, Babi, Khepry, and Imhotep, as defined by El-Maarry 
et al. (2015). They are located in the northern hemisphere, at latitudes between 0 and +60 
degrees and longitudes between 60 and 150 degrees. Because the data set covers sev- 
eral morphological regions, we expect significant variations in roughness values. Facets 
shared by the two datasets cover mainly the northern part of Imhotep and a small fraction 
of Aten. 


VIRTIS radiance at 4.95,1m A VIRTIS radiance at 4.95,1m A 


2014-09-12T23:28 2014-09-15T02:51 
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Figure 5.3: Radiance measurements from VIRTIS projected onto the shape model on 12 
Sep 2014 (left) and 15 Sep. 2014 (right). 


In order to compare data and model, we first convert the modelled surface tempera- 
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ture into radiance, /(A, T), by combining reflected solar spectral radiance and emitted IR 
radiation as follows 
I(A,T) 2 aF(A) * eB(A, T) (5.7) 


where a and e are the constant geometric albedo and emissivity (respectively 0.05 and 
0.95), F the solar flux modulated by a reflectance function (the Lommel-Seeliger Law), 
and B the black body radiation at wavelength A and temperature T. We neglected the 
directionality of epsilon, which is affected by spatially unresolved surface roughness. 

For a given facet on the shape model, the solar spectral radiance is given by Planck's 
law, heliocentric distance, and a photometric model. We chose the Lommel-Seeliger 
(Hapke 2012) correction, which depends only on incidence and emission angles. This 
equation is written as 





2 Rsa) 2h 1 
Fa) - x(? $ | ge (5.8) 


Ho tH 
where h,c,k are Planck’s constant, the speed of light, and the Boltzmann constant, 
respectively. The value Rs. is the sun radius, r, the heliocentric distance. The values uo 
and u are the cosines of incidence and emission angle. 
Finally, the emitted IR radiation is the Planck function at T and A, is 


h 
Th JD ensfmx-] 


2 
Bde ge (5.9) 
oA exr —] 

We applied our model to all facets covered by the VIRTIS data set, except those in 
shadows for which a more complex photometric model should be considered, accounting 
for multiple scattering. Shadows were calculated with a custom ray tracing engine! widely 
used to simulate Rosetta images, and calibrated against observations. 


5.4.5 Importance of roughness 


In addition to TI and emissivity, it has been noted by many authors that roughness at 
spatial scales unresolved by the instrument plays an important role when determining 
the surface temperature of an airless body with a low TI, at many length scales (Kührt 
et al. 1992, Groussin et al. 2013, Keihm et al. 2012, Davidsson et al. 2015). First, the 
micro-topography (on centimetre scale) of the surface leads to re-radiation and multiple 
scattering of the incoming solar and thermal radiation, which can locally enhance or de- 
crease the amount of energy received by forming heat traps or sinks. The temperature 
therefore varies on subpixel scales and we consequently observe a thermal radiation that 
cannot simply be represented by a single Planck function. One example is the case of 
craters, in which the shadowed wall can be heated up by radiation from the illuminated 
wall. Second, the illuminated fraction of the surface seen by the detector is influenced 
by the roughness, and so the energy received is a function of the relative position of the 
observer and the energy source. 

Both issues limit the accuracy of temperature retrieval by remote sensing and must be 
accounted for. It should be noted that simulating accurately the roughness at the subpixel 
resolution of the OSIRIS (Keller et al. 2007, Optical, Spectroscopic, and Infrared Remote 





'shapeViewer: www.comet-toolbox.com 
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Imaging System) camera over a large surface is computationally far too expensive to be 
achieved in most cases. In practice, the surface temperature is first approximated with a 
classical Lambert surface. In a second step, the roughness is modelled by multiplying the 
Lambert temperature with a correction factor that is a function of wavelength, incidence 
angle, emission angle, projected phase angle, and a statistical description of the subpixel 
topography. 

The subpixel topography can be represented in many ways, including fractal surface, 
coherent noise, and craters. Davidsson et al. (2015) showed that, for non-directional ter- 
rains (i.e. with facets orientated in all directions), it does not matter much in the final result 
which description is used. We adopted the approach proposed by Kührt et al. (1992), rep- 
resenting topographic variations with a set of spherically shaped depressions. The model 
is scale independent, and therefore our only free parameters are the area fraction covered 
by the craters and opening angle of the spherical segments. The opening angle refers to 
the centre of the corresponding sphere so an angle of 180° creates a half-spherical depres- 
sion. 

In this study, we want to assess whether the roughness of surface of 67P can be es- 
timated reliably from the combination of infrared measurements by the Rosetta VIRTIS 
instrument and a thermal model of the surface, described in sections 5.3 and 5.4.2, re- 
spectively. To do so, we assume that after accounting for the TI, any further discrepancy 
between model and data can be corrected by a factor that results solely from roughness 
effects. In this approach, we first derive a best TI for the observed area and then allow 
roughness to vary spatially in our model. We know all the geometric parameters for any 
facet of the shape model at any time, so we can use the roughness model to relate this cor- 
rection factor to the area covered by the spherical segment depressions. This roughness 
estimation is then compared to OSIRIS images of the same region, as they have typically 
much higher spatial resolution than VIRTIS data. We finally assess how the spatial dis- 
tribution of our retrieved thermal roughness (micrometre - millimetre scale) compares to 
the visual roughness (centimetre scale) in the images. 

We therefore have three unknowns: TI and two roughness parameters (crater density 
and opening angle). Ideally combining three inputs would provide good constraints on 
the parameters: for instance the radiance in at least three wavelength channels for a given 
facet, or radiance measurements at the same wavelength but acquired with different ge- 
ometries. However, since this is not always available in our dataset, we choose a less 
constrained approach, as outlined above. 


5.5 Results 


5.5.1 MIRO results 


Fig. 5.4, 5.5, 5.6, 5.7, and 5.8 show the SBT as functions of the electrical skin depth 
for each of the thermal inertia values used in the thermal model, in the range 5 to 320 
JK^!m^?s 9^, for the five dates given in Table 5.2. As stated earlier, we constrain the TI 
when SBT intersects MBT for electrical skin depths less than ~3.0 m and ~1.0 m for the 
millimetre and submillimetre cases, respectively. 

For example, in Fig. 5.4 in the millimetre channel, the SBT with TI = 5 JK ! m ?s ?? 
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Figure 5.4: Simulated brightness temperatures (SBT; coloured lines), as functions of the 
electrical skin depth, for various values of TI in the range 5 - 320 JK ! m ?s ?? (ST units). 
The black dotted lines represent the MBT from the millimetre and submillimetre receivers 
for September 1, 2014 at 23:48. The millimetre SBT curves with TI values less than 
10 JK-! m^s"9? intersect the millimetre MBT line; and the submillimetre SBT curves 
with TI values less than 80 JK^! m^?s ^9? intersect the submillimetre MBT line. 


intersects the MBT but TI 2 10 JK^! m ?s ?? does not. We hence derive that the TI in the 
millimetre channel for September 1 is «10 JK ^! m ?s ?^, as given in Table 5.2. 

We choose not to further constrain the TI value by interpolating between the curves. 
Linear interpolation seems to be insufficient, so instead of devising a complex interpola- 
tion scheme, we quote the upper and lower bounds (where possible) from the discrete TI 
values, which initially went into the thermal model. 

The MBT of the different observations performed with the millimetre receiver are con- 
sistent with a TI less than 80 JK ! m ?s ?^, while the value for September 1, 2014 implies 
that the TI should be less than 10 JK!m ?s 9^, The results from the submillimetre re- 
ceiver on September 1 and 2, 2014 are also consistent with a TI less than 80 JK^! m^?s 92, 
but the other dates imply higher values of 160 and 320 JK^! m ?s 9? are also possible for 
their observations. The results are compiled in the Table 5.2. 


5.5.2  VIRTIS results 


Disentangling TI from roughness effects is a complex problem for which we do not always 
have enough constraints. To reiterate, we attempt to fit first the TI as it relates directly to 
the physical properties of the material. We then add the surface roughness as a second 
parameter effectively describing which subpixel area contributes to the flux that reaches 
the detector. 
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Figure 5.5: Same as Fig. 5.4 but for September 2, 2014. The results are given in Table 
5:2. 
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Figure 5.6: Same as Fig. 5.4 but for September 12, 2014. The results are given in Table 
3.2. 
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MIRO/mm Observation, 13-09-2014-09:33 
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Figure 5.7: 
5:2. 
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Same as Fig. 5.4 but for September 13, 2014. The results are given in Table 
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Same as Fig. 5.4 but for September 15, 2014. The results are given in Table 
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Date Time Thermal inertia (JK 'm ?s ??) Region 
from mm observation from smm observation 
01/09/2014 23:48 «10 «80 Seth 
02/09/2014 04:23 «80 «40 Aten 
12/09/2014 23:28 «80 5-320 Ash 
13/09/2014 09:33 <40 <160 Ash 
15/09/2014 02:51 <80 >20 Aten 





Table 5.2: Range of possible local TI derived from the comparison between the SBT and 
the MBT for each MIRO observation. 


Therefore, the first step in our inversion of roughness parameters is to determine which 
thermal inertia brings our model the closest to the measured radiance. We calculated the 
radiance for the seven values of TI investigated by the thermal model described above 
(TI = 5, 10, 20, 40, 80, 160, 320 JK-!m^?s-??), We found that the modelled flux is 
almost always below the measured flux (see Fig. 5.9). As we increase TI, the average flux 
decreases. This is because higher TI means that facets must be illuminated for a longer 
time before reaching their maximum temperature. Our simulations show that a good fit is 
achieved for a TI larger than 80 JK^!m^?s ??. We therefore use this TI value as a lower 
limit in our analysis. Higher values lead to an equally good fit, suggesting that TI is not 
the main reason for the dispersion we observe in Fig. 5.9, although the fit is slightly worse 
for TI values of 160 and 320 JK^! m^?s 9? in the case of September 15, 4.75 um. 
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Figure 5.9: Left: Example of modelled radiance minus observed radiance at 4.95 um 
on September 15, 2014 for TI-5 JK^!m?s ??. Increasing the TI makes the plot more 
symmetrical around zero. Right: Mean value of model minus observation as a function 
of TI for two epochs and two wavelengths. We observe that increasing the TI beyond 80 
JK^!m^?s^?? does not improve the fit, and slightly worsens it for the case September 15, 
4.75 um. We therefore use this TI value as a lower limit in our analysis. 


Having selected TI=80 JK~'m~’s~°° as our reference TI, we then postulate that the 
remaining discrepancy between model and observation is caused by the surface rough- 
ness. Following the approach of Kiihrt et al. (1992), we define the corrected flux for a 
rough terrain as Z, = c, X I, with c, a function of incidence and emission angles defined 
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by the shape model facet normals, projected phase, crater density, crater opening angle, 
and wavelength. The projected phase, or azimuth, is the angle between the observer di- 
rection and solar direction vectors projected onto the facet surface. In this approach, the 
only free parameters are the crater density and the crater opening angles. Our test cases 
show that the influence of opening angle is less important than that of the crater density 
for the wavelengths considered. Therefore, we fix the opening angle to a value of 180 
(which leads to maximal c,) and find the crater density that minimizes the quantity Z, — Im: 
the difference between the calculated radiance, /,, and the VIRTIS measured radiance, 
In. The crater density is allowed to vary between 0 and iG the maximum density for 
packing equal circles on a two-dimensional plane. This approach effectively gives us a 
lower bound for the crater density on each facet in the region of interest. 

We note that in principle, one should find the best parameters by varying both TI and 
roughness in parallel when getting the minimum least-squares fit. However, this is a large 
parameter space and running the complete inversion goes beyond the scope of this work. 
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Figure 5.10: Difference between modelled rough radiances and VIRTIS measurements 

for September 12, 2014 (left) and September 15, 2014 (right). Although we still observe 

some scattering, the dispersion is far smaller than when not considering roughness, as in 

Fig. 5.9. 


Fig. 5.10 shows an example of the corrected radiances for two different dates. The 
most remarkable effect is the strongly reduced scattering of our data points with respect 
to Fig. 5.9, indicating a very good agreement with the observations. The residual dis- 
crepancies can be explained by the fact that we have used a relatively coarse shape model 
(124938 facets, 20 m resolution) and the same value of TI for all facets. A finer spatial 
resolution for both parameters would certainly improve the fits, as would a more ad- 
vanced roughness model and spatially varying emissivity. In addition, it should be noted 
that VIRTIS measurements are not instantaneous, but require up to 20 minutes in this 
work, during which the viewing geometry changes as the comet rotates and the spacecraft 
moves. This also affects our inversion. 

Figures 5.11 and 5.12 illustrate our results; the retrieved crater density is projected on 
the shape model. One can see that the distribution of densities is not random, but rather 
organized in regions whose distribution resembles the morphological regions derived from 
optical data alone. 
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Figure 5.11: Left: Regional map of the comet from El-Maarry et al. (2015) with morpho- 
logical units. Right: Roughness parameter, i.e. the density of the half-spherical depres- 
sions projected onto the shape model. The density value is represented by the colour bar 
on the right. 





Figure 5.12: Half-spherical depression density overlaid on OSIRIS images, 3D projected 
onto the shape model with the same colour scheme as in Fig. 5.11. The images are 
centred on approximately 130° longitude and 10° latitude. One can see that areas with 
similar morphologies share similar roughness. Visually smooth terrains (at the scale of 
OSIRIS data, i.e. 75 cm/px) can display high thermal roughness, indicating that the 
spatial scale of this parameter is far below the resolution of our images. 


We find a similar crater density for regions with a similar morphology, although this 
retrieved parameter may appear sometimes counter-intuitive. For instance, when looking 
at morphology alone, the Imhotep plain appears to be separated in two areas with different 
textures. The northern part is rougher with many metre-size boulders, while the southern 
area is very smooth and has topographic variations that are basically undetectable at the 
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best OSIRIS resolution (<50 cm/px). Interestingly, our roughness inversion also discrim- 
inates between these two subregions, thus helping to validate our method, but the area 
with highest crater density is the subregion that appeared to be the smoothest in OSIRIS 
images. 

Although this may seem surprising at first, it is actually explained by the fact that 
the roughness that plays a role for IR measurements is simply at a much lower scale 
than topographic variations. Taking this into account, it seems reasonable to postulate 
that consolidated terrain is likely to show uniform textures on a centimetre scale, while a 
visually smooth terrain may be covered in pebbles and dust aggregates that form a very 
rough surface at small scales. 

This result is consistent with several observations from the Rosetta lander Philae. Im- 
ages acquired by the ROLIS camera during the descent to the Agilkia region on November 
12, 2014 show that a region that appeared as a smooth dusty plain in OSIRIS images is 
actually covered by centimetre-sized pebbles when observed at higher resolution (Mottola 
et al. 2015). This region has a very granular appearance as seen by ROLIS (Rosetta Lan- 
der Imaging System) and may be a result of physical processes occurring on the comet 
such as airfall of small material or talus deposits. The same may be true for the Imhotep 
dust plain, as revealed by our thermal analysis. On the other hand, CIVA (Comet Infrared 
and Visible Analyser) images obtained at Abydos after landing show that the consoli- 
dated terrain, although very rough on a metre-size scale, display flat and smooth facets at 
centimetre to decimetre scale (Bibring et al. 2015). Flat and smooth morphologies could 
be a result of a variety of process, including the disaggregation of boulders by complete 
erosion to create rounded structures, or the fragmentation of hard consolidated material 
to reveal flat surfaces (Mottola et al. 2015). Further study is required to link the measured 
roughness to these morphological processes. 


5.6 Discussion and conclusions 


5.6.1 Thermal inertia 


All of the results obtained from the MIRO millimetre observations imply that the TI has 
a low value, «80 JK^!m^?s ??. The result from September 1 even implies that it should 
be «10 JK^! m ?s ?? in the Seth region. For this date, the submillimetre channel gives a 
higher upper bound of 80 JK ^! m ?s ?^, similar to other results. 

Constraining the TI from the subsurface temperature observations is sometimes dif- 
ficult, however. For the measurements on September 12, 2014 and September 15, 2014, 
the submillimetre observation cannot constrain the TIso well (5-320 JK-!m^?s^?? and 
220 JK ! m ?s 9^, respectively) based on the assumed dielectric model. The results from 
the millimetre observations on these two dates are better constrained however, giving an 
upper bound of 80 JK^! m?s-95. 

In comparison, VIRTIS work implies a best-fitting value of 80 JK !m ?s 9? across 
the observed Aten, Babi, Khepry, and Imhotep regions. The VIRTIS value was calculated 
over a much larger area than the MIRO values as the viewing area of VIRTIS is larger 
than the MIRO beam area. However, this does not seem to make a difference. With the 
VIRTIS data restricted to a subset containing only the facets viewed by the MIRO beam, 
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the VIRTIS result remains the same; a similar thermal inertia value is predicted. 
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Figure 5.13: Summary of measured TIs (units: JK !m ?s 9?) from a variety of sources. 
Measurements in black represent previous published works and coloured bars represent 
this work. From left to right, these measurements represent the TI from: (a) Spitzer obser- 
vations (Lamy et al. 2008); (b) a thermophysical model applied to light curve observations 
(Lowry et al. 2012); (c) MIRO observations by Gulkis et al. (2015); (d) MIRO observa- 
tions by Schloerb et al. (2015) (e) MIRO observations by Choukroun et al. (2015), with 
TI from the millimetre channel in black and TI from the submillimetre channel in grey; (f) 
MUPUS observations, with average represented in white (Spohn et al. 2015); (g) VIRTIS 
observations (Leyrat et al. 2016); (h) MIRO observations from this work in Seth region, 
with TI from the millimetre channel in dark green and TI from the submillimetre channel 
in light green (horizontal lines represent upper or lower bounds); (1) MIRO observations 
from this work in Aten region; (j) MIRO observations from this work in Ash region; (k) 
MIRO observations from this work in Ash region; (I) MIRO observations from this work 
in Aten region; and (m) VIRTIS observations from this work; the best-fitting value is 
indicated in white. 


In Figure 5.13, we compare our measurements for the TI from this work with pre- 
vious observations. Our measurements for the TI from MIRO are consistent with past 
MIRO measurements of the TI (Gulkis et al. 2015, Schloerb et al. 2015, Choukroun et al. 
2015) and there is also some overlap with the lander observation (Spohn et al. 2015). The 
measured VIRTIS TI is in good agreement with the MUPUS value, which was also de- 
termined from an infrared instrument: MUPUS-TM operates at 5 - 25 um and VIRTIS 
observes at 2-5 um. In addition, further VIRTIS calculations have yielded a TI in the 
range 25 - 170 JK^! m ?s ?? (Leyrat et al. 2016). Neglecting spatial and depth variations, 
a bulk TI value in the range 30-50 JK! m ?s ?? would fit all of the measurements from 
Rosetta instruments, except for the MIRO millimetre measurement made here on Septem- 
ber 1. This value is similar to the early TI estimates from Lamy et al. (2008) and Lowry 
et al. (2012), which are much lower than most of the Rosetta instrument measurements. 
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Taken together, the Rosetta measurements are consistent with a low TI 
(<170 JK"'m~’s~°>), and there is notable overlap in the error bars of these results. 


5.6.2. Roughness 


In the MIRO analysis, roughness is considered on the scale of the horizontal shape model 
resolution because the angle between each facet normal of the shape model and the di- 
rection to the spacecraft is taken into account when computing the SBT. We considered 
adding subfacet roughness into the MIRO analysis but this would be computationally 
intensive. In addition, the uncertainty on the TI is larger than the expected uncertainty 
caused by subfacet roughness because of the valid range of values considered for the 
electrical skin depth. 

The VIRTIS analysis showed that the roughness could be considered on subpixel 
scales that are below the resolution of the OSIRIS images. While work remains to be 
done to define the surface roughness properly, our approach shows that one can use IR 
radiance measurements and a simplified roughness model to separate the surface into re- 
gions of different properties. This study raises the prospect of consistently mapping the 
full surface at various epochs, in particular in areas like Imhotep in which significant 
changes have been detected at large and small scales over the course of the Rosetta mis- 
sion (Groussin et al. 2015). We conclude that the roughness below the resolution of the 
images plays a key role in interpreting thermal infrared observations, and that one cannot 
use the visible images to estimate the surface roughness at this scale, since the small- 
scale (IR wavelength) roughness is not well correlated with the roughness on the large 
scale (visible wavelength). 

This work demonstrates the importance of comparative studies of results from dif- 
ferent instruments. By combining observations from the MIRO, VIRTIS, and OSIRIS 
instruments on Rosetta, we have gained new insights into the small-scale roughness of 
the surface and found more evidence for a low TI of the comet nucleus. 
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6.1 Introduction 


Ever since Whipple (1950), we have understood comets as small solar system bodies 
made from volatile ices and refractory material. As they fly close to the sun, comets heat 
up as a result of the increase in received flux and sublimate their volatile components, 
forming a tenuous gas atmosphere around the nucleus called the coma. 

One of the initial goals of the MIRO instrument was to use the acquired spectral lines 
to study the evolution and development of the inner coma (Gulkis et al. 2007). The inner 
coma is unobservable from ground-based telescopes, and in situ instruments can only 
probe this region if they are in it. MIRO offers a unique opportunity to observe this region 
around comet 67P/Churyumov-Gerasimenko remotely, and throughout the duration of the 
Rosetta mission. 

Spectral line inversion has long been applied to Earth for meteorological purposes to 
obtain temperature profiles in the atmosphere (for example, King (1952)) and now more 
recently to other solar system bodies (Hartogh et al. (2010a) for Mars; Lellouch et al. 
(2015) for Uranus and Neptune; Thelen et al. (2018), Rengel et al. (2014), Bauduin et al. 
(2018) for Titan). As a result of the Rosetta mission, this can now be applied to a comet 
for the first time and the parameters of interest in the coma are the gas kinetic temperature, 
the expansion velocity (always given as the velocity of the gas towards the spacecraft) and 
the number density of water molecules. The H5!6O line is optically thick in the MIRO 
observations whilst the H5^O line is optically thin (Gulkis et al. 2015) and the different 
opacities of each line give sensitivity to different altitudes of the coma. 

Currently, the inner coma region is typically modelled with parametrised functions 
(Lee et al. 2015, Marshall et al. 2017, Hansen et al. 2016) but the actual behaviour may be 
far more irregular. Inverting the observed spectral lines can reveal the turbulent dynamics 
of cometary activity which parametrised profiles cannot capture. Activity is well ob- 
served but unfortunately not well understood. For instance, it is not always known which 
molecular species is driving the outgassing at a given time, how much of the surface is 
active, how the activity distribution changes, or even how the non-volatile components of 
the comet are lifted from the surface. By thoroughly characterising the inner coma, we 
can hope to address some of these open questions as the behaviour at the surface may be 
reflected in the dynamics of the coma. 

The aim of this chapter is to introduce the application of inversion techniques to the 
MIRO spectra which can then be used in future analysis of the data. In the next section, 
I will describe the method of iterative regularised inversion and in Section 6.3.1, I will 
apply this to a synthetic case study to test the accuracy of the profile retrieval. This will 
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then be applied to MIRO data in Section 6.3.2. 


6.2 Method 


In this section, I introduce the method which will be used for the inversion of the MIRO 
data. The method uses a regularised iterative inversion process called Optimal Estimation 
(Rodgers 2000) which starts from an a priori profile to constrain the solution. A single set 
of line of sight profiles for gas temperature, expansion velocity and number density are 
retrieved from a single pair of H56O and H5'5O observations. 


6.2.1 Creating the a priori profiles 


The first step is to create the a priori profiles from which the inversion will be constrained. 
Since there is very little prior knowledge available for the conditions of the inner coma, 
this is done by utilising the atmospheric information in the nadir water spectral lines: 


e The continuum temperature, Te, measured by MIRO can act as proxy for the near 
surface temperature in the coma. In reality this is a sub surface temperature which 
may be different to the temperature of the gas in the lowest coma altitudes 


e The temperature of the flat saturated H;!°O line core, T ,, gives an estimate for the 
coma temperature closest to the spacecraft altitudes 


e The Doppler shift of the H5'*O line, Vexp, Can be used in the a priori as the projected 
velocity in the coma around altitudes close to the spacecraft. 


e The line area ratio between the H,!°O and H5'*O lines can determine an estimate 
for the column density, N, as described in Marshall et al. (2017) which can then be 
used to set the integrated area under the density profile. 


This information from the spectral line helps to create a priori profiles for the tem- 
perature, velocity and density profiles. The temperature and velocity a priori profiles are 
spline functions which are schematically represented in Figure 6.1. Both spline functions 
have 5 knots given in Table 6.1. 

The first two knots are at 5.5 m and 100 m and take temperature and velocity values 
T,, and v,,, respectively. These are chosen within the limits T, - 10 K X T4, < Te + 10 K 
and 0 km/s 2 v,, 2 —0.4 km/s. T, is measured from the MIRO spectra as mentioned 
above. One of the major constraints when creating the a priori profiles is the lack of 
information about the velocity close to the nucleus. There is not much information about 
the velocity flow in this region in the MIRO spectra so our a priori profiles take a range of 
near-nucleus velocities from 0-0.4 km/s. This range of values covers possible values for 


the most probable thermal speed, v, = NS = y922T, where R = 8.3 J mol! K^! is the 
gas constant, M = 0.018 kg mol ! is the molar mass of water, and T; is assumed to be in 
the range 0-200 K. 

The fifth knot is at the altitude of the spacecraft H;;. and the temperature and velocity 
in this region, T',,. and v,,., respectively, are chosen within the limits T, - l5 K € Tj, X 
T; * 15 K and vo, + 0.1 km/s 2 vj; 2 voy — 0.1 km/s. 
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Figure 6.1: Schematic depiction of the creation of the a priori profiles for temperature and 
velocity. The a priori are spline functions with five knots (red circles). The co-ordinates of 
each knot are randomly generated from the intervals given in Table 6.1. The dashed lines 
show how different a priori profiles can be constructed where T, and v,, are the temper- 
ature and velocity close to the nucleus, 7. and v; are the values close to the spacecraft, 
and 7,;; and v, are intermediate values. The temperature and velocity values are esti- 
mated from the measured parameters of the MIRO spectra (continuum temperature, 7, 
the saturated H5'9O line core, T,, and the Doppler shift of the H2'%O line, v..)). 
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Knot | Altitude (m) | Temperature (K) | Velocity (km/s) 
1 5.5 Tgr Ver 
2 100 ie Ver 
3 500-800 Forte verte 
4 1000-5600 Te Vs/c 
5 He. Ts LO 




















Table 6.1: Co-ordinates for the knots of the a priori spline functions of the temperature 
and velocity profiles. T,, is the ground temperature close to the nucleus randomly chosen 
from the interval T. — 10 K € T,, < Te + 10 K, where T, is the measured continuum 
temperature. T/c is the spacecraft temperature randomly chosen from the interval 7, — 
15 K < Tye < Ts + 15 K, where T, is the measured temperature of the saturated H,!°0 
line core. v,, is the expansion velocity at the nucleus randomly chosen from the interval 
O km/s 2 vg, 2 —0.4 km/s. vs,. is the spacecraft expansion velocity randomly chosen 
from the interval vo; + 0.1 km/s > Vsje = Vexp — 0.1 km/s, where Vex, iS the (negative) 
measured Doppler shift expansion velocity in the direction of the spacecraft. Hsje is the 
spacecraft altitude. The altitudes for knots 3 and 4 are randomly chosen from the given 


intervals. 


The fourth knot has the same values for temperature and velocity as the fifth knot, 7',;. 
and v,,., but this time, the altitude of this point is randomly chosen between 1000-5600 m. 


Finally, the third knots for temperature and velocity are randomly chosen to be in 


the interval 500-800 m and they take values of 7,,;; and v,;; where Tj; = a and 
rt ^s. v 
Vmid — "s T / . 


The co-ordinates of the knots are randomly chosen within their respective intervals for 
two main reasons. The first is that the measured atmospheric parameters from the MIRO 
data (T.e, T, and v,,;) are only proxies for the atmospheric values that they represent. In 
addition, these parameters have their own error bars. 

Secondly, by setting the knots within these limits, a series of a priori profiles can 
be generated. The inversion can then be performed multiple times, each time starting 
from a different a priori profile with different coordinates for the knots. In this way, we 
can quantify the dependence of the solution on the a priori parameters and hence assess 
the error on the retrieval. The a priori temperature and velocity profiles are randomly 
created by choosing T,, 7’; and v,,, from within their uncertainty ranges and by randomly 
selecting the third and fourth knot altitudes. The variety in a priori starting profiles will 
be useful later when we want to quantify the error on our retrieval solutions. 


The shape of the density profile is estimated from DSMC modelling of the inner coma 
region (Marschall et al. 2016), behaving with an approximate 4 dependence above a few 
hundred metres and having a reasonably constant density below this altitude. The density 
profile is scaled so that the integrated column density under the curve is equal to the 
column density, N, as estimated from the line area ratios of the H2!6O and H5!*O lines. 
The density profile is only scaled in the inversion process so only one scaling parameter 
is retrieved to describe the behaviour of the density. The grey dotted lines in Figure 6.3 
give an example of the a priori profiles generated in this way. Again, there is uncertainty 
on the measurement N, so the a priori density profile has a randomly chosen integrated 
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column density within the error bars of N - a factor of about 2-3 (Marshall et al. 2017). 

Each profile is interpolated to have 31 grid points and thus 30 cells. The first two 
grid points are at altitudes of 5.5m and 100m and after that, the remaining 29 points are 
distributed logarithmically until the altitude of the spacecraft. Using a purely logarithmic 
grid gives too many altitude points very close to the nucleus (within 100m) where the 
sensitivity in the measurements is not so high. 


6.2.2 Creating K,S, and S, 


The next step involves creating the weighting functions or kernel, K, the covariance matrix 
Sa, and the measurement error covariance matrix, S e. 

K is made by perturbing the elements in the a priori profiles. For the velocity and 
temperature profile, a small perturbation (1 K in temperature and 20 m/s in velocity) is 
made at each altitude grid point and the density profile is perturbed by a single scaling 
factor (0.1). In this way, 61 perturbed profiles are made and each one is run through a 
radiative transfer equation in local thermal equilibrium to simulate the water spectral lines 
seen by MIRO. Each simulated pair of spectral lines is subtracted from the simulated lines 
of the unperturbed a priori profile to see how each perturbation has altered the spectral 
lines. The difference between the perturbed and unperturbed lines are stacked in a matrix. 
This is the kernel. It effectively shows how each perturbation produces changes in the 
spectral lines. This will allow the inversion process to know at which altitude changes 
can be made in the profiles to minimise the difference between the a priori spectra and the 
observed spectra. The LTE radiative transfer method used here is the same as that used 
in Lee et al. (2015) and Marshall et al. (2017), using an approximated escape probability 
approach (Bockelée-Morvan 1987, Litvak and Kuiper 1982). 

S, takes the form of a square diagonal matrix with the $.[i,i] values assigned to be 
the error on each point of the spectra, e;, and all other values are set to zero. The matrix 
looks like: 


e; 0 0 0 
0 e 0 .. 0 

Sg=]0) O vez a 0 (6.1) 
00 0 .. e 


and it has dimensions of m x m, where m is the spectral resolution of the water lines. 

Sa is also a square matrix but it has dimensions of p x p where p is the number of 
perturbations through the three atmospheres. For single beam inversion, p — 61. The 
diagonals of the covariance matrix effectively describe how constrained the inversion is 
to the a priori. Off diagonal elements allow for consecutive altitude elements to be corre- 
lated. S, is not a diagonal matrix, and we instead use a Gaussian distribution centred on 
the diagonal element. The elements of S, are found using the equations: 








l |j- ky ! 
Sali, k] = T2 exp] -" 7 l; forj.k <31 (6.2) 
cl 
| 02 
SER Y. ew[-2 m jor 31 € j k «61 (6.3) 
cl 
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S „[61,61] = D},, (6.4) 


for each set of profiles. Ty:a, Vsrq and D,,. are the standard deviations of the temper- 
ature and velocity, or how much each variable is allowed to deviate from the a priori. 
Ta and Vy are the respective correlation lengths or how many altitude grid points are 
correlated together, and j and k are elements of S,. For the inversions performed here, 
T sta =3 K, Viu = 0.02 km/s, D a = 0.1, T — 5 and Va =5; 


6.2.3 Optimal estimation inversion 


In the third step, I use a regularised optimal inversion technique from Rodgers (2000), 
incorporating K, S, and S, as described above. The equation takes the form: 


Xii — Xi * G[y - F(xi) * K(Gi — xa)] (6.5) 
G=(KS,K’+S,)" (6.6) 


where i denotes the iteration, G is called the gain matrix, y is the observed spectra, 
F(xj) is the simulated spectra from the x; profiles, x, is the a priori profile, and K? is the 
transposed kernel. In Equation 6.5, the difference between observed and simulated spectra 
is minimised, and regularised by the smoothness of the profiles and the measurement error 
on the observed spectra. In the first iteration, x; = x, but after that x; changes with each 
iteration. The process now starts again, taking the x;,, retrieved solution and creating 
a new kernel by perturbing the elements in these new profiles. Typically, the inversion 
converges on an answer in 5-10 steps. 

It is possible to calculate an error covariance matrix from our inputs but because there 
are some unknown constraints in our inversion (the model parameter error, the forward 
model error), we prefer to quantify the error by examining the non-uniqueness of starting 
from each different a priori. By using the error on the measured parameters mentioned in 
Section 6.2.1, we create a family of a priori profiles and start the inversion from each one. 
Each a priori gives a different final answer and the variation in the retrieved answers gives 
us the error on each parameter. 

For single beam inversion, a single set of temperature, velocity and density profiles are 
retrieved along a single line of sight. The a priori parameters are taken from the observed 
spectra. x;, x, and x;,, have dimensions of 1 x 61, y and F(x;) have dimensions of 1 x p 
where p — 260 for the spectral line resolution used here for the MIRO observations. K isa 
61x260 matrix and S , and S , are square matrices with 61 and 260 elements, respectively. 


6.3 Application 


In order to apply the inversion method to observed MIRO data, it is necessary to first 
check the performance of the parameter retrieval on a synthetic case study. In this way, 
the accuracy of the retrieved solutions can be quantified against prior determined true 
profiles. The geometry used for the synthetic case study and the MIRO case study is 
shown in Figure 6.2 and this is taken from April 9, 2016 at 23:48 when MIRO was 29.9 km 
away from the nucleus. The MIRO beam is shown in red and the temperature, velocity 
and density will be retrieved along this line of sight. 
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Figure 6.2: Left: Geometry of the MIRO observation on April 9, 2016 at 23:48. The red 
line shows the length of the MIRO beam from the spacecraft to the nucleus (29.9 km) and 
the yellow lines shows the direction of the sun. Right: zoomed in view of the geometry. 
MIRO is observing a part of the relatively flat Imhotep region which is well illuminated 
at this time. 


6.3.1 Synthetic case study 


For the synthetic case study, a coma atmosphere is created with profiles for temperature 
velocity and density. These are shown as black lines in the top panel of Figure 6.3. From 
these ‘true’ profiles, synthetic MIRO measurements are made for the H5'^O and H5;'*O 
lines. These are shown in black in the top panel of Figure 6.4. Noise is added to the signal 
at a level comparable to that of an observation made with 12-minute integration time, as 
will be used for the MIRO data in the next section. The aim now is to retrieve the true 
synthetic profiles from the synthetic water line spectra using the inverse method described 
above. 

The first step is to estimate the coma parameters Te, Ts, Ve», and N from the synthetic 
MIRO spectra. These are found to be: 26 K, 185 K, -0.76 km/s and 4.12x10P? cm. 
These values are then used to create the 40 a priori profiles shown in grey dashed lines 
in the top panel of Figure 6.3. From these a priori profiles, the inversion begins and 
Equation 6.5 is solved. The retrieved solutions are shown in the yellow dashed lines in 
the top panel of Figure 6.3 and the difference between the retrieved profiles and the true 
profiles are shown in the middle panel of Figure 6.3, again in yellow dotted lines. The 
average of the retrieved solutions is then shown in green for the top and middle panels. 

Qualitatively, the temperature profile and the velocity above 1 km is retrieved very 
well, and the average solution follows the true profiles, finding the inversion in both pro- 
files. The solution for the density profile is also good, being reasonably close to the true 
profile. The difference between the average retrieved solution and the true profiles con- 
verges well around the zeroth line, shown in grey in the middle panels of Figure 6.3. The 
retrieved spectra are shown as coloured lines in Figure 6.4 and the residuals, the difference 
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Figure 6.3: Top: Retrieved profiles for gas temperature, expansion velocity and water 
number density from the synthetic MIRO observations shown in Figure 6.4 using the 
viewing geometry of Figure 6.2. The black lines are the ‘true’ synthetic profiles, the grey 
lines show the a priori profiles, the yellow lines show the retrieved answers from each a 
priori and the green line is an average of all of these answers. Middle: The difference 
between the true and retrieved profiles for each inversion is shown in the yellow dotted 
lines. The green line is the average of all of these. Bottom: The absolute difference 
between the retrieved and true profiles are shown in yellow. The average of all of these 
profiles is shown in green and this quantity is called the bias. The standard deviation of 
the bias is called the variability and it is shown in grey. 
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Figure 6.4: The synthetic MIRO spectra of H56O and H5!*O (black lines) and the re- 
trieved spectra (coloured lines) along with the residuals (bottom panel). 
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between the retrieved and measured spectra are all within 8 K, with the biggest differences 
correlating to the noisier parts of the spectra. 


The bottom panels of Figure 6.3 give the absolute difference between the average 
retrieved solution and the true solution, and from this we will quantify the quality of the 
retrieval by defining the errors. The yellow dotted lines again show the absolute difference 
between the true and retrieved profiles for each retrieval. The green line then gives the 
average of all of these. We call this the bias. For the temperature profile, only a small 
region at around 200 m has a bias greater than 10 K and in fact above 1250 m from the 
nucleus, the bias is less than 8 K. The bias in the velocity profile is less than 30 m/s from 
the true profile above 1250 m and often much better, particularly at around 3000 m. The 
bias below 1000 m is not so good however, but this is to be expected. The MIRO water 
line spectra are not particularly sensitive to the velocity close to the nucleus. There is a 
small bias in the density profile, underestimating the density by a factor of about 8%. 


We also define the variability which is the standard deviation of the bias. It measures 
the variability in all of the retrievals. It is shown in grey in the bottom panel of Figure 
6.3. Generally the 1 c variability is close in value to the bias except in temperature less 
than 1000 m from the surface where the variability in the retrieved profiles decreases in 
comparison to the bias. We will use the bias to define the error on our retrieval from the 
real MIRO spectra. 


6.3.2. MIRO case study 


Now that we have assessed the performance of the inversion method, it can now be ap- 
plied to real spectra measured by MIRO. Figure 6.5 shows the retrieved profiles for gas 
temperature, expansion velocity along the line of sight and water number density in the 
top panel, and the MIRO spectra underneath captured between April 9, 23:48 and April 
10, 00:00, 2016. From the MIRO spectra, the measured coma parameters (Te, Ts, Vexp and 
N) are: 41 K, 190 K, -0.66 km/s and 1.04x10'° cm ?. The grey a priori lines in the top 
panel of Figure 6.5 were made from these values. The retrieved solutions are shown in the 
yellow dashed lines and the average solution is given in green. The retrieved spectra are 
again shown as coloured lines in the bottom panel of the figure along with the residuals. 


The retrieved solutions are reasonably consistent with little variability in the temper- 
ature and density profiles or the velocity profile above ~1000 m. The temperature profile 
increases from 190+5 K close to the surface (5.5 m) to 196+9 K at 111 m before de- 
creasing to a minimum temperature of 50+3 K at 17200 m. The temperature then stays 
reasonably constant up to the altitude of the spacecraft. The integrated column density 
under the number density profile takes a value of (2.33+0.19)x10!° cm”. The velocity 
at 1000 m is -0.54+0.04 km/s and increases to a value of -0.70+0.01km/s at an altitude 
of 4200 m. There is then a decrease towards a velocity minima of -0.67+0.02 km/s at 
11500 m before finally increasing again to -0.70+0.03 km/s at the altitude of the space- 
craft. The error on each parameter is taken from the bias given in the bottom panel of 
Figure 6.3. 
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Top: Retrieved profiles for gas temperature, expansion velocity and water 


number density from the MIRO lines observed on April 9, 2016 at 23:48 with the viewing 
geometry shown in Figure 6.2. The grey lines show the a priori profiles, the yellow 
lines show the retrieved answers from each a priori and the green line is an average of 
all of these answers. Bottom: The MIRO spectra (black lines) and the retrieved spectra 
(coloured lines) along with the residuals. 
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6.4 Discussion 


In this section, we have shown how an iterative optimal estimation inversion method can 
be applied to MIRO data to retrieve coma profiles for gas temperature, expansion velocity 
and water number density. This method started from a priori profiles which were made 
using coma parameters estimated from the MIRO spectra, and through the regularised 
inversion method, I arrived at solutions for the profiles in the coma. 

For the MIRO observation made on April 9, 2016 at 23:48, the temperature can be 
seen to decrease steadily from approximately 190 K close to the nucleus to 50 K at the 
spacecraft. The velocity increased from 0.54 km/s at 1000 m to a maxima of 0.70 km/s at 
before decreasing again to a minima of 0.67 km/s at 11500 m. The velocity then appeared 
to increase once more with a value of 0.70 km/s close to the spacecraft. The water had a 
column density of 2.33x10P cm? for this observation. 

The inversion method was tested with a synthetic case study to assess the errors on the 
parameters. The density profile has an error of about 8% at all altitudes as this profile has 
only one scaling parameter. Above 1250 m from the nucleus, the errors on temperature 
and velocity are quite small, less than 8 K and less than 30 m/s. Below this level, the error 
on the temperature can increase to 14 K mainly owing to the fact that the temperature 
profile is slightly overestimated at the altitudes where the temperature decreases rapidly. 
This can happen when the density is slightly underestimated as in this example. Below 
1000 m, the velocity profile is poorly retrieved. The MIRO spectra has little sensitivity 
here so the velocity cannot be found close to the nucleus. 

The application of inverse methods will enable the dynamics of the inner coma region 
to be properly studied for the first time and allow structures like inversions to be discov- 
ered in the temperature and velocity profiles. By establishing the gas behaviour, it will be 
possible to characterise the temporal and spatial variations in the inner coma and study 
the activity in this unexplored region. Assessing the activity will help to answer some of 
the open questions about the physics of cometary outgassing. 

Deriving coma profiles will help with comet modelling as it will enable the coma to 
be described from a data driven perspective for the first time, rather than with assumed 
parametrised profiles. In addition, the open questions about activity - distribution, driving 
species, dust lifting mechanism - can potentially be addressed by characterising the gas 
flow around the nucleus. Finally, physical aspects like photolytic heating or a Knudsen 
layer (Maali et al. 2016) could potentially be identified from the profiles which would 
enhance our understanding of the physical behaviour of outgassing in the coma. 

These questions will be addressed in the future and not only from the nadir observa- 
tions described here. It will also be possible to invert the spectra from limb observations. 
In addition, utilising the geometry of observations could potentially improve the accuracy 
on the retrieval, as is done in medical tomography. These will require new methods to be 
developed and tested before they can be applied to the MIRO data. 
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The science objectives for the MIRO instrument were outlined in Section 1.6 and these 
aims included: determining the abundances of major volatile species; characterising prop- 
erties of the nucleus; and exploring the behaviour of the inner coma close to the surface 
(within 100 km). Through this work, I have shown how MIRO is uniquely placed to pro- 
vide constraints on the nucleus and the coma of comet 67P/Churyumov-Gerasimenko and 
achieve the goals set out before the beginning of the Rosetta mission. Here, the results 
are briefly summarised and put into a wider context. 


Water abundance 


In Chapter 3, I used the observed water spectral lines to measure the water production 
rate pre- and post-perihelion from the line area ratios of the H5'^O and H5!*O lines (Mar- 
shall et al. 2017). The abundance of sublimating water increased with decreasing he- 
liocentric distance, as expected for an icy body, and reached a measured maximum of 
(1.4220.51)x10?? molec/s or (426 + 153) kg/s on August 29, 2015. For comparison, the 
measured production rates 231 days before perihelion and 245 days after perihelion were, 
(3.24+1.79)x10°° molec/s, and (1.72+0.54)x10?° molec/s, respectively. The peak pro- 
duction from linear fitting of the data suggests that the peak production occurred 34410 
days after perihelion. 

During the perihelion passage, the | MIRO results indicate that 
67P/Churyumov-Gerasimenko lost (2.4 + 1.1)x10° kg of water, or (1.2 + 0.6)x10!° kg of 
total material, assuming a dust-to-gas ratio, y = 4 (Rotundi et al. 2015). Taking the mass 
of the comet to be 1.0x 10P kg (Sierks et al. 2015), this gave a mass loss of 0.12 + 0.06 %. 
With a small mass loss estimated for this apparition, it is tempting to assume that much of 
the nucleus interior is unprocessed during its lifetime, and the content beneath the surface 
may still be pristine. However, many aspects of the comet are still unknown. The dynam- 
ical history of 67P is not well determined more than 200 years ago and it is possible that 
it had close perihelion passages with the sun before (Maquet 2015), thus processing the 
material in the comet. In addition, the collisional history is also not well known and nei- 
ther are the mechanisms which cause outbursts, possibly driven by volatile material in the 
deep interior. Comets are complex bodies and more must be done to properly understand 
their activity before conclusions can be drawn regarding solar system history. 

The ROSINA instrument also provided an opportunity to measure the water produc- 
tion rate from Rosetta, using the onboard pressure sensor and mass spectrometer to de- 
termine the abundance of water around the spacecraft (Hansen et al. 2016). With a direct 
simulation Monte Carlo model, they find higher production rates than reported here. One 
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reason for this may be the different acquisition techniques of the two instruments. MIRO 
is a remote sensing instrument and so observes material along its line of sight whereas 
ROSINA measures the material in the immediate vicinity. Differences between the re- 
sults from both instruments may come from how both methods treat the behaviour of the 
gas between the nucleus and spacecraft, and raises the possibility of non-r? dependence on 
the density, possibly due to extended sources. The discrepancies between measurements 
from the two instruments needs to be explained in future work. 

In addition, the spatial resolution of MIRO enabled the water production to be limited 
to a determined region on the surface. The surface was found to be quite inhomogeneous 
in terms of outgassing, with Wosret, Neith and Bes being the regions where the highest 
measured production rates originated from. Notably, these are some of the southernmost 
regions and are particularly well illuminated around perihelion. MIRO did not observe 
high production rates from regions on the northern hemisphere, except for a couple of 
measurements in the Seth and Ash regions. 

I also found the gradients for the change in the water production to have values of 
-3.8 and -4.3, pre- and post-perihelion, respectively. So far though, no thermal modelling 
which is consistent with the observed MIRO water densities has been able to reproduce 
the detailed variation in Q(r;) and, as discussed in Chapter 4, the meaning of these slope 
values is not entirely clear. Crucially, the gradients of the production rate curves are also 
a function of heliocentric distance and a single value cannot describe the outgassing be- 
haviour over large heliocentric distances. It is nearly impossible to disentangle the effects 
of nucleus shape, rotation axis orientation and activity distribution from the water pro- 
duction rate curve if these parameters are unknown, and when these three parameters are 
known as for 67P, the question now turns to how activity works to produce the observed 
outgassing. The inability of a simple sublimation model to reproduce features of the 67P 
water production curve suggest that there is inhomogeneous activity distributed on the 
surface but the cause of the inhomogeneity is still unknown. It could be explained by in- 
homogeneous ice coverage across the surface or due to dust coverage which has the effect 
of suppressing volatile sublimation. 


Nucleus properties 


Using the continuum channels of MIRO, I estimated the thermal inertia of the nucleus, 
as described in Chapter 5 (Marshall et al. 2018). In the context of all other Rosetta mea- 
surements, it seems likely that the surface has a low average value for the thermal inertia 
with 50 JK~'m~s~°° fitting within the error bars of most of the observations from MIRO, 
VIRTIS and MUPUS. Looking more widely, the derived values for other comets for the 
thermal inertias can also be low, with values less than 50 JK^!m^?s 9? for unresolved 
nuclei (Groussin et al. 2018). Results from 67P (Marshall et al. 2018), 9P/Tempel and 
103P/Hartley (Davidsson et al. 2013, Groussin et al. 2013) imply that higher thermal in- 
ertias may also be possible. This could be explained by variations in composition across 
the surface or by the fact that thermal inertia is a function of temperature, which varies 
across the surface based on the solar illumination. Comets appear then to have lower 
thermal inertias than many other solar system bodies, including near-Earth and Main Belt 
asteroids (Dell’Oro et al. 2007, Tanga et al. 2009). Looking to other solar system bod- 
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ies, the Saturnian satellites (Howett et al. 2010), Pluto (Lellouch et al. 2017) and Trans- 
Neptunian Objects (Lellouch et al. 2013) appear to have very low thermal inertias too, 
less than 25 JK ^! m?s ?? (Groussin et al. 2018). 

As a result of low values for the thermal inertia, the temperature in the comet surface 
decreases quickly with depth. Depending upon the modelling used, after a couple of 
metres the temperature in the nucleus can be in the range 30-100 K (Marshall et al. 2018, 
Hu et al. 2017). The erosion rate during an orbit is on the order of decimetres or less 
(Gortsas et al. 2011, Keller et al. 2015, Groussin et al. 2018). These two findings imply 
that much of the nucleus could be unprocessed since formation if the interior parts of 
the comet remain very cold, and that samples from beneath the surface could be pristine 
remnants of the early solar system. 

By combining this analysis with the VIRTIS data, it was also possible to map the 
small-scale thermal roughness on the shape model. 


The inner coma 


The application of inversion methods as described in Chapter 6 enables the retrieval of 
gas temperature, expansion velocity and water number density along the MIRO line of 
sight from spectral line observations. In this chapter an example inversion was performed 
for observations made on April 9, 2016. Starting from a priori profiles derived from coma 
parameters in the MIRO spectra, I arrived at profiles which show the steady decrease in 
temperature and the increase of expansion velocity from 1000 m above the surface to the 
spacecraft. I was also able to find the column density of water in the coma. 

In the future, the widespread application of inversion methods will enable the inner 
coma to be probed and characterised in unprecedented detail and give new insights into 
the outgassing behaviour of comets. Through the analysis of the MIRO data, the temporal 
and spatial evolution of the coma activity can be studied and hence help to further our 
understanding of how activity works. Activity is generally not well understood - what 
drives it, which species dominate and how the material is lifted. Additionally, the retrieved 
profiles can help to establish the flow field around the nucleus. It may then be possible 
to assess the physics of the coma like a potential Knudsen layer close to the surface. The 
profiles could also reveal any photolytic heating or extended sources of water in the coma. 

It is also possible that certain viewing geometries may be able to help improve the 
accuracy of the retrieval approach. Several observations made over a short period time 
which view the same volume of the coma from different perspectives could be inverted 
together to retrieve a single set of parameter profiles. This would be similar to tomograph- 
ical methods employed in medical studies. Additionally, consecutive measured could be 
correlated together, thereby further regularising the problem. 


Concluding remarks 


MIRO provided an opportunity to probe the activity of comet 
67P/Churyumov-Gerasimenko as part of the Rosetta mission and it has already revealed 
inhomogeneities in the spatial and temporal outgassing from the surface and in the coma. 
Much of the MIRO dataset has not been investigated yet, with most results so far coming 
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from the early part of the mission. As a result, the MIRO dataset will be a valuable re- 
source to use in any attempts to properly understand the physical behaviour of cometary 
activity. Hopefully there are still many exciting discoveries to come from the comet. 
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Appendix 


A Derivation of radiative transfer 
equation 


This is the derivation of the equation used for microwave remote sensing given in Equation 
from the standard definition of the radiative transfer equation (Equation 2.8). We start with 
this equation 


dI,(s) 
ds 
and divide both sides by a,(s) 


] dh) js) 





= JvCS) Em a, s)MyCs) (A.1) 











= -I, A2 
as) ds a,(s) iu oy 
This can then be expressed as 
I, 
N 5-19) (A3) 
dr 


using the definition of the source function and the fact that dr= a,(s)ds. Allthe terms 
can then be multiplied by exp(r) 


dl,(s) 








exp(T) dr t L(s)exp(r) 2 S,(s) exp(r) (A.4) 
which can be simplified using the Product Rule 
di, d 
dpt MOOS Oe (A.5) 
dr dr 
to get 
d 
a ^09 exp(t) = S,(s) exp(t) (A.6) 


Writing /,(5) exp(r) and S,(s) exp(t) as 1, and $,, respectively, this now becomes 


di, . 
es; A 
dr ( ) 


This differential equation can be solved using 


I, TO 
a di, = [ $, dr (A.8) 
ho 0 
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where the right side can simply be redefined in terms of S,(s) and the left hand side 
becomes 


A 


I, 
ii di, = AR = 1, - 1,0 = I,expfr} - 1,9 exp{0} (A.9) 
I, 


$ I, v0 
‚0 


which then leads to the form of the radiative transfer equation 


LE T Lo exp(-T,(So)) + J- S, exp(-T,($)) a,(s) ds (A.10) 
0 


that is commonly used in microwave remote sensing. 
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